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I.  INTRODUCTION 

The  possibility  that  a  large  asteroid  may  impact  the  Earth  delivering  an  energy  in 
excess  of  30  MT  of  TNT  is  a  very  real  threat.  One  only  has  to  look  at  the  Moon  in  the 
night  sky  to  gain  an  appreciation  for  the  magnitude  and  probability  of  such  an  impact. 
How  the  human  race  deals  with  the  threat  of  such  an  impact  is  of  significant  interest  to  the 
planet  Earth  as  a  whole. 

As  discussed  by  Rather  et.  al.  (1992),  as  early  as  1705,  when  Edmund  Halley 
wrote  A  Synopsis  of  the  Astronomy  of  Comets,  there  has  been  speculation  that  an 
extraterrestrial  object  might  impact  the  Earth,  the  possibility  of  such  impacts  was  not 
perceived  as  a  threat  to  life  on  Earth  until  the  late  1940's  when  the  Moon's  craters  were 
fully  understood  to  be  the  result  of  impact  events,  not  volcanism,  and  that  the  Earth  is 
subject  to  the  same  impact  hazard.  Rather  et.  al  (1992)  also  mention  that  the  magnitude 
of  the  threat  was  not  fully  realized  until  the  late  1970s  to  early  1980s.    The  perceived 
threat  began  to  stand  upon  a  solid  foundation  with  the  publication  by  Alvarez  et.  al. 
(1980)  of  a  theory  on  the  extinction  of  the  dinosaurs  due  to  the  impact  of  a  large  asteroid 
or  comet  65  million  years  ago.  The  theory  put  forth  by  Alverez  et.  al.  suggests  that  an 
impact  by  a  10  km  diameter  asteroid  at  the  Chicxulub  site  off  of  the  Yucatan  peninsula 
indirectly  caused  the  extinction  of  60%  of  all  life  on  the  Earth,  including  the  dinosaurs. 

Since  the  awareness  of  the  possibility  of  an  asteroid  impacting  the  Earth  began 
developing  in  the  mid  20th  century,  there  have  been  several  "near  misses"  recorded. 
Perhaps  the  most  spectacular  near  miss  was  a  large  fireball  created  by  an  object  racing 
across  the  daytime  sky  in  a  northerly  direction  that  entered  the  Earth's  atmosphere  above 
Utah  in  the  United  States  and  exited  above  Alberta,  Canada  on  the  10th  of  August,  1972. 
The  object  observed  was  determined  to  be  an  asteroid  upwards  of  30  meters  in  diameter 
as  reported  in  Sky  and  Telescope  Magazine  (1972).  Had  this  asteroid's  trajectory  been 
ever  so  slightly  different,  mankind  would  have  had  its  first  opportunity  to  observe  a  large 
object  impact  the  Earth.  Such  an  impacting  object  would  carve  out  a  crater  200  to  500 
meters  in  diameter.  Since  then,  several  other  asteroids  and  comets  have  been  detected 
passing  by  the  Earth  at  distances  less  than  a  few  hundred  thousand  kilometers.  Apollo 
Asteroid  1989FC,  referenced  by  the  AIAA  Space  Systems  Technical  Committee;  Asteroid 
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1991BA,  reference  by  Scotti  et.  al  (1991);  and  Asteroid  1996JA1  referenced  by  Jaroff 
(1 996)  are  three  asteroids  discovered  recently  passing  very  close  to  Earth.    On  an 
astronomical  scale,  these  close  approaches  are  essentially  impacts. 

In  April  of  1990  the  American  Institute  of  Aeronautics  and  Astronautics  issued  the 
position  paper  entitled  Dealing  with  the  Threat  of  an  Asteroid  Striking  the  Earth  that 
briefly  described  the  implications  of  an  asteroid  impact.  The  AIAA  found  that  "Earth- 
orbit-crossing  asteroids  clearly  present  a  danger  to  the  Earth  and  its  inhabitants."  It  was 
recommended  that  a  "systematic  and  open  program"  for  detection  of  Earth  crossing 
asteroids  be  established  as  well  as  a  study  to  "define  systems  which  can  deflect  or  destroy, 
or  significantly  alter  the  orbits  of,  asteroids  predicted  to  impact  the  Earth."  A  few  search 
programs  existed  prior  to  the  position  paper  and  a  few  others  have  begun  subsequently. 

As  a  result  of  the  awareness  of  the  possibility  of  an  asteroid  or  comet  impacting 
the  Earth,  the  National  Aeronautics  and  Space  Administration  has  sponsored  two 
workshops  to  study  the  fundamentals  of  the  impact  and  impact  mitigation  problem.  The 
first  workshop  is  summarized  in  the  report  entitled  "The  Spaceguard  Survey:  Report  of 
the  NASA-International  Near-Earth-Object  Detection  Workshop."  The  second  workshop 
is  summarized  in  the  report  entitled  "Near-Earth-Object  Interception  Workshop."  The 
concentration  of  the  workshops  have  been  related  to  assessing  the  magnitude  of  the  threat, 
impact  effects  and  hazards  to  the  Earth,  as  well  as  the  political  implications  of  developing 
an  impact  mitigation  capability.  Several  books  have  been  published  on  the  matter,  both 
technical  and  non-technical,  and  more  than  one  Hollywood  movie  has  been  based  on  the 
subject.  Two  spacecraft  exploration  missions,  Near  Earth  Asteroid  Rendevouz  (NEAR) 
and  Clementine,  have  included  intercepts  of  asteroids  as  a  major  part  of  their  mission  to 
study  the  nature  of  asteroids.  However,  little  "astrodynamical  analysis"  has  been 
performed  on  the  feasibility  of  the  impact  mitigation  problem. 

The  astrodynamical  analysis  of  feasibility  is  where  this  thesis  is  intended  to  fit  into 
the  larger  problem.  Presented  is  an  analysis  of  the  impact  and  mitigation  problem  based 
on  a  two-dimensional  and  two-body  analysis.  It  is  intended  to  be  a  first  order 
approximation  for  the  solution  of  a  larger  problem.  However,  it  is  a  more  rigorous 
treatment  of  the  astrodynamics  of  the  hazard  than  previous  published  analyses,  such  as 


Ahrens  and  Harris  (1994).  Included  in  the  following  analysis  is  the  periodic  nature  of  the 
problem  as  well  as  the  near  term  effects.  The  analysis  is  not  concerned  with  object 
detection  or  orbit  prediction,  but  instead  centers  on  how  impulses  applied  to  an  asteroid  at 
various  points  on  the  asteroid's  orbit  affect  the  outcome  when  there  is  a  presumption  of 
collision  otherwise.  Mission  design  for  mitigation  is  not  a  goal  of  the  following  analysis; 
however,  the  analysis  tool  presented  may  be  utilized  in  determining  a  first  order  estimate 
for  optimizing  the  time  and  position  of  asteroid  intercept  for  impact  mitigation. 

Presented  first  is  an  astronomical  development  of  the  existence  of  an  asteroid 
impact  problem  from  the  origins  of  the  solar  system  to  the  record  of  past  impacts  on  the 
Earth  followed  by  a  brief  description  of  two-body  orbits.  The  problem  is  then  presented 
with  governing  assumptions  and  a  method  of  solution.  The  solution  method  is  then 
assessed  and  validated  followed  by  an  analysis  of  an  asteroid  impact  scenario  with  a 
discussion  of  the  results.  The  analysis  method  is  then  applied  to  the  asteroid  Toutatis 
which  will  make  several  close  approaches  to  the  Earth  over  the  next  decade. 


II.         SOLAR  SYSTEM  OBJECT  ORIGINS 

A.  SOLAR  SYSTEM  FORMATION 

Current  theories  of  formation  of  the  Solar  System  stem  from  the  post  initial 
expansion  universe  environment  of  gas,  dust,  radiation,  and  magnetic  fields  in  a 
nonuniform  distribution.    Bouyed  by  the  interplay  of  gravitational,  magnetic,  and  pressure 
forces  local  mass  concentrations  began  to  form.  This  interplay  of  forces  on  the  non- 
uniform environment  sets  up  initial  angular  momentum  conditions  for  very  large  three- 
dimensional  structures.  These  large  rotating  structures  began  to  coalesce  into  accretion 
disks  around  central  masses.  These  central  masses  eventually  reached  critical  mass  for 
fusion  to  occur  and  stars  emerged.  We  call  our  local  star  the  Sun. 

B.  PLANETARY  FORMATION 

In  a  similar  fashion  to  the  stellar  formation,  the  accretion  disk  around  the  Sun 
provided  the  environment  for  smaller  mass  concentrations  and  accretion  disks  to  form  thus 
producing  small  scale  structures  called  planetismals  composed  of  the  basic  chemical 
elements  in  varying  quantities.  Gravitational  attraction  and  relative  motion  of  these 
planetismals  caused  them  to  collide  with  one  another  and  cohesive  forces  enabled  some  to 
remain  attached  forming  larger  structures.  After  enough  of  these  interactions  occurred, 
the  planets  began  to  form.  Unlike  for  the  stellar  conditions,  the  planets  do  not  possess  the 
critical  mass  to  initiate  and  sustain  a  fusion  reaction.  This  allows  for  large  scale  assembly 
of  solid,  aqueous  and  gaseous  structures  to  take  place  in  quantities  and  composition 
proportional  the  relative  percentages  of  chemical  elements  present  in  the  planetismals. 
These  structures  are  more  familiar  on  the  Earth  as  the  crust,  the  oceans  and  the 
atmosphere.  Analysis  of  the  elements  present  from  the  impact  delivery  mechanism  and 
current  known  compositions  of  asteroids  and  comets  suggests  that  this  was  the  mechanism 
of  organic  and  non-organic  material  delivery  that  formed  the  Earth  as  proposed  by  Chyba 
et.  al.  (1994). 

C.  ASTEROID  AND  COMET  FORMATION 

However,  not  all  of  the  accretion  disk  surrounding  the  Sun  coalesced  into  either 
the  Sun,  the  Planets  or  their  satellites.  The  remainder  of  the  planetismals  continue  to  be 
dispersed  throughout  the  Solar  System  in  the  form  of  asteroids  and  comets.  The  asteroids 


being  located  primarily  within  the  inner  Solar  System,  and  the  comets  existing  in  the  outer 
Solar  System  and  beyond. 

The  asteroids  are  mainly  concentrated  in  the  asteroid  belt  located  between  the 
orbits  of  Mars  and  Jupiter.  The  remainder  of  the  asteroids  are  dispersed  throughout  the 
solar  system  in  elliptical  orbits  lying  mainly  inside  the  orbit  of  Jupiter.  There  are  several 
theories  of  how  the  asteroid  belt  came  to  exist.  These  ideas  include  the  destruction  of  a 
planet,  or  the  inability  of  a  planet  to  form,  due  to  the  tidal  forces  of  Jupiter.  There  are 
several  theories  accounting  for  the  formation  of  the  asteroids  not  located  in  the  asteroid 
belt.  These  ideas  range  from  planetismals  that  have  never  collided  and  attached 
themselves  to  another  planetary  body,  to  cast  off  remnants  of  massive  collisions  of 
planetary  bodies  with  very  large  planetismals,  to  objects  from  the  asteroid  belt  that  may 
have  been  perturbed  by  a  passing  object  into  a  smaller  orbit.  All  of  the  ideas  have  some 
amount  of  merit  that  give  them  validity. 

Comets  are  believed  to  originate  from  the  Oort  Cloud  of  cometary  material 
orbiting  the  Sun  at  a  distance  of  some  50,000 AU.  From  this  cloud,  comets  are  believed  to 
be  injected  into  the  solar  system  by  orbital  perturbations  due  to  the  gravitational  field  of 
passing  stars.  Once  injected  into  the  solar  system,  gravitational  encounters  with  the  Sun 
and  Jupiter  may  further  perturb  the  comets'  orbit  and  either  "capture"  the  comet,  so  it 
remains  within  the  solar  system,  or  "assist"  the  comet  on  its  way  to  interstellar  space. 
Additionally,  there  is  believed  to  exist  a  band  of  icy  objects  that  extends  from  the  orbit  of 
Neptune  at  30  AU  out  to  as  much  as  100  AU.  These  objects  are  said  to  be  located  in  the 
Kuiper  Belt,  so  named  after  Gerard  Kuiper,  who  first  proposed  their  existence  in  1 95 1 . 
The  objects  that  form  the  Kuiper  belt  range  in  size  and  orbital  characteristics  to  the  extent 
that  it  is  believed  that  Pluto  may  actually  be  a  member  of  this  group  as  discussed  by  Jewitt 
and  Luu  (1995). 
D.         ASTEROIDS 

1.  Location 

The  asteroids  are  of  particular  interest  as  potential  impact  hazards  in  that  they  have 
a  greater  mass  density  than  comets  and  are  more  likely  to  reach  the  surface  of  the  Earth  in 
a  given  impact  scenario.  The  majority  of  the  asteroids  are  located  in  the  asteroid  belt 


between  the  orbits  of  Mars  and  Jupiter.  Very  few  asteroids  have  been  detected  beyond  the 
orbit  of  Jupiter.  This  lack  of  detection  may  only  be  an  effect  of  the  limited  capability  of 
current  detection  sensors.  However,  a  significant  number  of  asteroids  in  orbits  smaller 
than  that  of  the  typical  asteroid  belt  object  have  been  identified.  These  are  the  objects  of 
primary  concern  for  the  problem  of  mitigation. 

2.  Quantities 

Literally  thousands  of  asteroids  have  been  observed  and  identified  orbiting  the  Sun. 
Of  those,  more  than  300  are  considered  near  Earth  asteroids  (NEA's)  and  pose  a  threat  as 
a  potential  impacting  object.  Of  greater  concern  are  the  subset  of  NEA's  dubbed  Earth 
crossing  asteroids  (ECA's)  which  are  currently  about  200  in  number.  The  orbits  of  the 
ECA's  are  such  that  they  allow  for  the  possibility  of  impact  with  the  Earth  at  some  future 
date.  These  ECA's  range  in  size  from  10  km  down  to  0.1  km,  which  is  currently  the  limit 
of  detection  of  asteroids  by  Earth  based  sensors.  Estimates  for  ECA's  near  and  above  10 
km  in  diameter  indicate  that  all  of  the  asteroids  have  been  identified  and  that  for  the  ECA's 
of  1  km  in  diameter  or  less,  the  identified  asteroids  represent  only  about  10%  of  those  that 
are  believed  to  exist,  as  stated  by  Grieve  and  Shoemaker  (1994). 

3.  Classification 

Classification  of  the  ECA's  are  determined  with  respect  to  the  Earth's  orbital 
extrema.  The  Earth's  perihelion  and  aphelion  distances  are  0.9833  and  1.0167  AU 
respectively.  The  orbits  of  the  ECA's  all  have  perihelia  less  than  the  aphelion  of  the  Earth 
orbit  and  aphelia  greater  than  the  perihelion  of  the  Earth  orbit.  The  ECA's  have  been 
subdivided  into  three  classes,  the  Atens,  Apollos  and  Amors,  based  on  their  orbital 
characteristics  as  discussed  by  Rabinowitz  et.  al.  (1994).  Table  1  summarizes  the 
distinction  between  the  three  classes.  Figures  1,  2  and  3  depict  typical  orbits  of  each  of 
the  three  classes. 


Class  Name 

Semi-Major  Axis 
a(AU) 

Perihelion  Distance 

q(AU) 

Aphelion  Distance 
Q(AU) 

Aten 

<1 

- 

>  0.9833 

Apollo 

>1 

<  1.0167 

- 

Amor 

- 

1.0167  <q<  1.3 

- 

Table  1 .  ECA  Classes 

2062  Aten,  Semi-major  axis  =  0.9666  AU,  Eccentricity  =  0  1826 


Figure  1 .  Typical  Aten  Orbit 


4179  Toutatis,  Semi-major  axis  =  2.5154  AU,  Eccentricity  =  0.6361 


Figure  2.  Typical  Apollo  Orbit 


433  Eros,  Semi-major  axis  =  1 .4582  AU,  Eccentricity  =  0.2229 


Figure  3.  Typical  Amor  Orbit 


4.         Physical  Properties 

The  physical  properties  of  the  asteroids  are  generally  that  of  "rocky",  irregularly 
shaped  spinning  objects.  Estimates  of  asteroid  densities  range  from  a  more  cometary 
density  of  2x1 03  kg  m"3  to  a  dense  asteroid  of  5x1 03  kg  m"3  with  a  mean  density  of  about 
3x1 03  kg  m3.  Asteroids  have  been  observed  that  are  composed  of  a  single  solid  mass  as 
well  as  multiple  mass  centers  either  physically  connected  or  gravitationally  bound  at  a 
contact  surface.  As  stated  by  Winters  (1996),  some  of  the  asteroids  may  actually  be 
aggregates  of  numerous  smaller  bodies  that  are  gravitationally  bound  together.  This 
hypotheses  is  further  supported  by  analysis  of  object  spin  motion.  It  appears  the  many 
asteroids  spin  at  angular  rates  sufficiently  slow  as  to  permit  gravitational  binding.  The 
case  of  the  comet  Shoemaker-Levy  9  supports  these  theories  in  that  it  was  gravitationally 
separated  by  tidal  forces  from  a  previous  passage  of  Jupiter  prior  to  its  1994  impact. 
E.         COMETS 

1.  Location 

Early  in  the  20th  century,  Jan  Oort  hypothesized  that  virtually  all  of  the  comets 
originate  in  a  cloud  of  cometary  matter  beginning  some  5  0,000 AU  from  the  Sun.  The  so 
called  Oort  cloud  is  assumed  to  be  a  uniformly  distributed  spherical  shell  that  surrounds 
the  solar  system  and  extends  to  approximately  half  of  the  distance  to  the  nearest  star, 
Alpha  Centauri.  From  this  cloud,  comets  are  injected  into  the  solar  system  by  orbital 
perturbations  of  passing  stars.  Once  inside  the  solar  system  the  comets  may  be  further 
perturbed  by  the  planets.  The  planetary  perturbations  may  "capture"  the  comet  in  the 
interior  of  the  solar  system  or  it  may  assist  the  comet  on  its  way  ejecting  it  from  the  solar 
system  for  all  time. 

2.  Quantities 

The  number  of  comets  in  the  Oort  cloud  is  believed  to  be  diminishing,  however  the 
remaining  number  is  believed  to  be  on  the  order  of  1012.  The  quantity  of  comets  that  exist 
within  the  solar  system  is  far  smaller,  only  100  or  fewer  have  been  identified. 

3.  Classification 

The  "captured"  comets  are  classified  as  periodic  comets  and  are  subdivided  into 
two  groups  as  summarized  by  Shoemaker  et.  al.  (1994).  Members  of  "Jupiter's  Family" 


10 


have  aphelia  close  to  Jupiter's  mean  orbital  distance  from  the  Sun  and  have  periods  of  less 
than  20  years.  Members  of  the  "Halley  family"  are  the  so  called  long  period  comets  and 
have  orbital  periods  greater  than  20  years  but  less  than  200  years. 

4.  Physical  Properties 

Comets  are  primarily  composed  of  "rock"  and  "ice."  As  such  they  have  been 
sometimes  called  "dirty  snowballs"  or  "icy  dirtballs"  depending  on  their  relative 
composition.  The  icy  material  is  believed  to  be  composed  of  water,  methane  or  ammonia. 
The  rocky  material  is  a  variety  of  carbonaceous  substances.  Their  mean  density  is  less 
than  that  of  the  asteroids  and  is  estimated  to  be  about  1000  to  2000  kg  m"3.  As  a  comet 
"burns"  off  its  icy  material  from  repeated  encounters  with  the  Sun,  the  nucleus  may  in  fact 
become  an  asteroid  of  a  somewhat  smaller  dimension. 
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III.       IMPACT  RECORD 

Occasionally  an  asteroid  or  a  comet's  orbit  is  such  that  it  actually  hits  another 
object  within  solar  system  such  as  the  event,  widely  celebrated  in  the  media,  of  comet 
Shoemaker-Levy  9  impacting  the  planet  Jupiter  on  the  16th  of  July,  1994.  The  very  same 
mechanism  that  caused  the  solar  system,  in  particular  the  Earth,  to  form  and  sustain  life 
now  threatens  the  very  same  life. 

A.  LUNAR  IMPACT  RECORD 

One  theory  of  the  origin  of  the  Moon,  as  mentioned  by  Chapman  and  Morrison 
(1994)  describes  it  forming  as  the  result  of  a  Mars  size  object  impacting  the  Earth  early  in 
its  existence  .  The  ejecta  from  the  event  is  believed  to  have  behaved  in  such  a  fashion  as 
to  remain  in  orbit  and  coalesce  into  what  is  now  the  Moon.  Further  evidence  of  such 
impacts  is  the  cratered  face  of  the  lunar  surface.  The  Moon  displays  literally  millions  of 
impact  sites.  Since  the  Moon  has  no  atmosphere  or  large  scale  geologic  processes  there  is 
no  erosion  of  the  impact  record,  unless  by  subsequent  impact  the  previous  impact 
structure  is  destroyed.  Thus,  the  Moon  serves  as  a  reminder  for  all  time  of  the  nature  of 
the  impact  hazard. 

B.  TERRESTRIAL  IMPACT  RECORD 

On  our  planet  Earth,  the  atmosphere,  oceans,  volcanism  and  plate  tectonics  tend  to 
erode  past  impact  sites.  As  discussed  by  Grieve  and  Shoemaker  (1994)  there  are  currently 
about  140  known  impact  craters  around  the  world.  The  locations  of  these  sites  are 
displayed  in  Figure  4. 
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Figure  4.  Known  Impact  Site  Locations  from  Grieve  and  Shoemaker  (1994) 

These  impact  craters  range  in  size  from  1 .5  km  in  up  to  200  km  in  diameter. 
Perhaps  the  best  example  of  a  classic  impact  crater  is  the  1 .5  km  diameter  Barringer  Crater 
located  in  Arizona,  see  Figure  5.  Barringer  Crater  is  believed  to  be  the  most  recent  impact 
site  on  Earth  having  formed  around  50,000  years  ago. 


Figure  5.  Barringer  Crater  from  Grieve  and  Shoemaker  (1994) 
However,  the  most  recent  impact  event  to  have  resulted  in  surface  damage  is 
believed  to  be  the  Tunguska  event  which  took  place  over  northern  Siberia  in  1908.  The 
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Tunguska  event  is  believed  to  be  the  result  of  a  60  m  diameter  comet  or  light  asteroid  that 
exploded  10  km  up  in  the  atmosphere  releasing  some  30  MT  of  energy  leveling  2500 
square  kilometers  of  forest.  The  Tunguska  event  is  characteristic  of  a  small  scale  impact 
occurrance.  On  the  large  scale  end  of  the  impact  spectrum  lies  the  K/T  impact  event  so 
called  by  its  occurrence  in  time  at  the  boundary  of  the  Cretaceous  and  Tertiary  periods.  It 
is  believed  that  an  enormous  asteroid  some  10  kilometers  in  diameter  created  the  200  km 
diameter  Chicxulub  impact  site  beneath  the  Gulf  of  Mexico  off  the  coast  of  the  Yucutan 
peninsula  and  is  responsible  for  the  extinction  of  60%  of  the  living  species  on  Earth  65 
million  years  ago. 

There  have  been  numerous  recent  close  calls  of  asteroids  impacting  the  Earth. 
Jaroff(1996)  describes  a  near  impact  as  recently  as  June  of  1996  where  an  object  roughly 
600m  in  diameter  passed  within  450,000  kilometers  (about  70  Earth  radii)  of  the  Earth. 

C.  OTHER  PLANETARY  IMPACTS 

Numerous  impact  sites  have  been  observed  by  planetary  exploration  spacecraft  that 
have  been  sent  to  Venus  and  Mars.  As  mentioned  above,  Comet  Shoemaker-Levy  9 
impacted  the  planet  Jupiter  in  July  of  1994. 

D.  IMPACT  SCALE 

The  energies  released  by  the  impact  of  asteroids  on  the  Earth  are  quite  large.  It  is 
conventional  to  express  the  impact  energy  in  terms  of  megatons  of  TNT  (MT),  where 
1MT  =  4.2x1 015  J.  With  the  assumption  of  a  mean  asteroid  density  of  3000  kg  m"3  a  size 
versus  impact  velocity  relation  may  be  made  for  a  given  impact  energy.  Figure  6  displays 
an  estimate  of  asteroid  mass  versus  impact  velocity  for  a  range  of  impact  energies  from 
tens  of  megatons  of  TNT  up  to  a  petaton  (1015)  of  TNT.  The  diameter  estimate  assumes 
an  effective  spherical  radius  corresponding  to  the  mass  for  the  same  impact  energy  and 
velocity. 
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Impact  Velocity  (km/s) 
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Figure  6.  Impact  Scale 

Estimates  of  impact  energy  and  terrestrial  devastation  have  lead  to  the 
classification  of  impacts  as  local,  regional,  and  global  events  as  described  in  Morrison  et. 
al.  (1994).  The  distinction  between  these  is  somewhat  blurred,  however  it  is  proposed 
that  local  events  correspond  to  impact  energies  in  the  vicinity  of  30  MT  or  less  which  will 
affect  approximately  0.001%  of  the  Earth's  surface  area  or  about  the  size  of  a  large 
metropolitan  area.  Regional  impact  events  are  those  that  release  energy  in  the  vicinity  of 
300  MT  to  3x1 04  MT  and  affect  about  0.1%  of  the  Earth's  surface  or  about  the  size  of  a 
large  state.  Global  events  are  considered  to  be  impacts  that  release  energies  near  and 
above  3x1 05  MT  affecting  approximately  10%  of  the  Earth's  surface  area  or  about  the  size 
of  a  large  country.  These  large  events  may  cause  such  disruption  of  the  ecological 
environment  by  particulate  injection  into  the  atmosphere  that  greater  that  25%  of  the 
Earth's  population  may  be  eliminated.  Figure  6  also  shows  several  representative  impact 
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scenarios  for  a  typical  impact  velocity  of  20  km  s" .  The  Tunguska  event  is  depicted  by  the 
'*'  symbol.  Tunguska  represents  a  nearly  classical  local  impact  event  of  an  object  some  60 
m  in  diameter.  The  recent  May  1996  near  miss  is  depicted  by  the  '+'  symbol  and 
represents  a  regional  event  for  an  object  some  600  m  in  diameter.  A  global  impact  is 
indicated  by  the  'o'  symbol  and  corresponds  to  an  asteroid  some  1 .5  km  in  diameter.  The 
'©'  symbol  represents  the  K/T  event  that  created  the  Chicxulub  impact  site  and 
corresponds  to  a  10  km  diameter  object.  It  is  noted  that  the  K/T  event  is  well  above  the 
threshold  for  global  catastrophe. 
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IV.       ORBITS 

A.         CONIC  SECTIONS 

In  a  simple  inverse  square  gravitational  field,  orbits  take  the  shape  of  conic 
sections.  For  an  object  that  lacks  sufficient  energy  to  escape  the  gravitational  attraction  of 
the  central  body,  the  resultant  orbit  will  take  the  shape  of  an  ellipse.  For  an  object  with 
sufficient  energy  to  escape  the  gravitational  attraction  of  the  central  body,  the  resulting 
orbit  will  take  the  shape  of  an  hyperbola.  The  transition  from  an  elliptical  to  hyperbolic 
orbit  occurs  along  an  "escape"  trajectory  shaped  as  a  parabola.  Straight  line  orbits  that 
lead  to  collision  of  the  orbiting  body  with  its  mass  center  are  special  cases  of  elliptical, 
parabolic,  and  hyperbolic  orbits. 


B. 


COORDINATE  FRAMES 


1.  Three-dimensional 

To  define  the  physical  problem  in  time  and  space  a  reference  frame  needs  to  be 
established.  For  the  general  three-dimensional  case  a  Cartesian  Sun-centered  inertial,  or 
Heliocentric,  coordinate  frame  is  chosen.  Figure  7  displays  the  orientation  of  the 
Heliocentric  coordinate  frame. 


first  day 
of  summer 


verna  I    equinox 
direction 


first   day 
of  winter 


( Seasons    are  for   Northern  Hemisphere  ) 


Figure  7.  Heliocentric  Coordinate  System  from  Bate  et.  al.  (1971) 
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The  direction  in  the  ecliptic  plane  from  the  Sun  to  the  First  point  of  Aries  serves  as 
the  primary  coordinate  direction  in  the  Heliocentric  frame.  The  ecliptic  normal  serves  as 
the  "vertical"  coordinate  and  is  defined  as  positive  in  the  northern  half  of  the  celestial 
sphere.  The  last  coordinate  direction  is  defined  by  taking  the  cross  product  of  the 
previous  two  coordinate  directions. 

2.  Two-dimensional 

To  simplify  matters  for  the  present  analysis,  a  two-dimensional  planar  perifocal 
coordinate  frame  is  chosen.  The  principal  axis  defining  a  two-dimensional  elliptical  orbit  is 
the  direction  toward  periapsis  from  the  primary  focus.  It  is  from  this  primary  axis  that  the 
true  anomaly  is  measured  in  a  counterclockwise  sense.  The  secondary  axis  is  normal  to 
the  primary  axis  in  a  right  handed  sense.  The  resultant  coordinate  system  is  shown  in 
Figure  8. 


*"  P ,  Perihelion 


Figure  8.  Perifocal  Coordinate  System 
C.         ORBITAL  ELEMENTS 

1.         Three-dimensional 

In  the  general  three-dimensional  case  six  orbital  elements  are  required  to  define  a 
particular  location  in  space  of  an  object  in  orbit.  Those  orbital  elements  are  the  semi- 
major  axis  (a),  eccentricity  (e),  inclination  to  the  ecliptic  (i),  longitude  of  ascending  node 
(Q),  argument  of  periapsis  (co)  and  the  time  of  periapsis  passage  (T).  See  Figure  9. 
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Figure  9.  Three-dimensional  Coordinates  from  Bate  et.  al.  (1971) 

2.  Two-dimensional 

In  the  perifocal  coordinate  system,  only  the  semimajor  axis  (a),  eccentricity  (e), 
and  true  anomaly  (v)  are  required  to  fix  a  position  on  an  orbit.  (Note  the  argument  of 
periapsis  is  defined  to  be  zero  in  this  coordinate  frame.).  Other  parameters  of  concern  are 
the  perihelion  distance  (q)  and  the  aphelion  distance  (Q).  Figure  10  displays  the  two- 
dimensional  perifocal  coordinate  system. 


»P,x 


Figure  10.  Perifocal  Coordinates 
D.         KEPLER'S  SECOND  LAW 

Motion  of  an  object  about  the  primary  focus  in  an  elliptical  orbit  is  governed  by 
Kepler's  Second  Law:  the  line  joining  the  planet  to  the  Sun  sweeps  out  equal  areas  in 
equal  times.  This  relation  increases  the  difficulty  in  the  problem  solution  in  that  there 
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exists  a  trancendental  relationship  between  time  and  position.  This  relationship  is  called 
Kepler's  Equation  and  takes  the  form: 

E-esinE 

t — — .  a) 

e  +  cos  v 

where  cos  E  = .  (2) 

1  +  ecosv 
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V.         PROBLEM  FORMULATION  AND  SOLUTION 

A.  STATEMENT 

Given  an  impending  asteroid  impact  with  the  Earth,  one  would  like  to  know  if 
there  is  an  optimum  point  on  the  asteroid's  orbit  where  an  impulse  may  be  applied  to  yield 
the  greatest  achievable  separation  distance  at  the  closest  point  of  approach  for  a  fixed 
impulse.  Additionally,  it  is  desirable  to  determine  if  is  there  an  optimum  direction 
associated  with  the  applied  impulse  that  further  increases  the  separation  distance. 

B.  ASSUMPTIONS 

In  order  to  proceed  with  an  approximate  solution  method,  simplifying  assumptions 
were  made.  The  first  assumption  was  that  of  two-body  motion.  That  is,  the  Sun  is  the 
mass  center  about  which  the  asteroid  and  the  Earth  orbit.  Furthermore,  the  asteroid  and 
Earth  do  not  gravitationally  interact  with  each  other.  It  was  also  assumed  that  there  were 
no  external  perturbing  effects  on  the  orbits  due  to  non-gravitational  forces  other  than  the 
applied  perturbing  impulse.  All  orbits  were  assumed  to  be  coplanar,  which  yields  a  two- 
dimensional  problem.  The  perturbing  impulse  is  assumed  to  occur  instantaneously.  The 
asteroid  is  assumed  to  be  one  of  the  near  Earth  objects  and  hence  in  an  elliptical  orbit 
around  the  Sun.  Hyperbolic  and  parabolic  orbits  were  not  considered  for  this  analysis. 
Finally,  it  was  assumed  that  the  Earth  is  in  a  perfectly  circular  orbit  at  1 AU. 

C.  TEMPORAL  CONSIDERATION 

In  determining  the  separation  of  a  NEO  from  the  Earth,  time  becomes  the 
dominant  factor  in  solving  the  problem.  The  relative  phase  of  each  of  the  orbiting  bodies 
determines  whether  the  bodies  will  collide.  Hence,  the  Earth-to-NEO  separation  distance 
is  the  quantity  of  concern.  Changing  the  orbital  elements  becomes  secondary  to  changing 
the  asteroid's  orbital  phase  with  respect  to  the  Earth. 

D.  METHOD  OF  SOLUTION 

Solution  to  the  above  problem  may  be  achieved  by  use  of  a  numerical  simulation 
scheme  where  the  orbital  equations  of  motion  are  integrated  from  some  initial  condition 
forward  in  time.  However,  such  a  simulation  is  time  consuming  (on  the  order  of  several 
minutes  per  impact  scenario)  and  therefore  limits  the  scope  of  any  analysis.  The 
assumptions  stated  above  allow  for  use  of  analytic  elliptical  orbit  equations.  Building  a 
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solution  based  on  these  analytic  two-body  equations  offers  a  greatly  reduced  time  for 
solution  (on  the  order  of  a  couple  seconds  per  impact  scenario)  and  therefore  a  greater 
scope  of  analysis  may  be  performed.  A  Mathworks  MATLAB  model  was  constructed  to 
numerically  execute  the  following  method.  A  numerical  integration  simulation  was  also 
developed  in  order  to  validate  the  analytic  method.  The  solution  description  below  is 
quite  general.  For  a  detailed  solution  description  see  Appendix  A.  The  MATLAB  script 
that  corresponds  to  the  solution  method  is  displayed  in  Appendix  B. 

1.  Geometry 

A  general  description  of  an  elliptical  orbit  intersecting  a  circular  orbit  is  used  in  the 
problem  solution.  This  description  suffices  for  all  planar  intercept  scenarios.  A  rotation  of 
the  perifocal  coordinates  may  be  required  to  bring  the  model  in  alignment  with  the 
particular  problem,  but  the  relative  geometry  remains  unchanged.  Figure  1 1  demonstrates 
this  equivalence. 


Figure  1 1 .  Equivalent  Impact  Scenarios 

To  uniquely  fix  an  intercept  scenario  with  the  above  assumptions,  only  the  impact 
true  anomaly  and  orbital  eccentricity  need  to  be  defined.  In  the  case  of  planar  orbits,  the 
perihelion  direction,  as  defined  by  original  asteroid  elliptical  orbit,  defines  the  principal  axis 
from  which  the  impact  true  anomaly  is  measured.  Implicit  in  the  above  impact  location 
description  is  the  assumption  of  an  Earth  orbital  radius  of  1  AU. 
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2.         Solution  Flowchart 

Figure  12  depicts  the  steps  in  the  problem  solution  method. 


Specify: 
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Figure  12.  Solution  Flowchart 
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3.  Given  Conditions 

For  a  given  impact  scenario,  it  is  assumed  the  impact  true  anomaly,  asteroid  orbital 
eccentricity  are  known  and  the  desired  impulse  time  prior  to  impact  is  specified.  It  should 
be  noted  here  that  specifying  impact  true  anomaly  and  orbital  eccentricity  fixes  the  semi- 
major  axis  and  hence,  the  orbital  period.  If  either  parameter  is  changed,  the  semi-major 
axis  changes.  That  is,  as  impact  true  anomaly  increases  from  0  to  n  the  semi-major  axis 
decreases,  thus  decreasing  the  orbital  period. 

4.  Unperturbed  Orbital  Elements 

From  the  given  conditions  and  the  known  Earth  orbit  the  perihelion  distance,  semi- 
major  axis,  and  orbital  period  for  the  unperturbed  orbit  may  be  determined  by: 

r^l  +  ecosvj 


l-e' 

and  P  =  — . 
n 


where  n  =   |— 


5.  Impact  Condition 

The  time  from  perihelion  of  the  impact  with  respect  to  the  unperturbed  orbit  may 
be  determined  by  Kepler's  Equation  as  given  in  Equations  (1)  and  (2). 

6.  Initial  Impulse  Condition 

The  time  from  perihelion  of  the  impulse  with  respect  to  the  unperturbed  orbit  may 
be  determined  by  Figure  13  and  the  equation: 

(Impulse  Time)Periheiion  =  (Impact  Time)penheiion  -  (Impulse  Time)impac,  Time- 
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(Impulse  Time)i„p.ctT,„ 


"-(Impulse  Tune)penb«i.on 

Figure  13.  Impulse  to  Impact  Time  Relation 


The  true  anomaly  at  impulse  may  be  determined  by  inverting  Kepler's  Equation 
such  that  true  anomaly  becomes  a  function  of  time.  This  approach  requires  an  iteration  on 
eccentric  anomaly. 

The  impulse  position,  f ,  and  velocity,  v ,  may  now  be  found  from: 

,(l-e=) 


r  = 


(3) 
(4) 


1  +  ecosv 
f  =  (r  cos  v)x  +  (r  sin  v)y , 

and  v  =  I— I  (-  sin  v)x  +  (e  +  cos  v)y  I , 

where  p  =  a(  1  -  e2  J . 

7.  Perturbation 

An  orbital  perturbation,  Av ,  is  applied  to  the  asteroid  at  the  impulse  position  r . 
This  yields  a  perturbed  orbital  velocity  v  +  Av . 

8.  Perturbed  Orbital  Elements 

From  the  impulse  position  f  and  perturbed  velocity  v  +  Av  the  perturbed  orbital 
elements  may  be  determined  from: 

h  =  f  x  (v  +  Av) , 
h2 


e  = 


P  = 


H 


H 


lv  +  Av|2--3-  r-{f  -(v  +  Av)}(v  +  Av) 
and  Equation  (4). 
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9.  Perturbed  Impulse  Condition 

Substituting  the  perturbed  orbital  elements  into  Equation  (3)  enables  the  true 
anomaly  at  impulse  to  be  determined  with  respect  to  the  perturbed  orbit.  The  time  from 
perihelion  of  impulse  with  respect  to  the  perturbed  orbit  may  be  found  from  Kepler's 
Equation. 

10.  Orbital  Positions  at  Impact  Time 

To  determine  the  separation  distance  of  the  perturbed  asteroid  from  the  Earth,  the 
positions,  and  times,  of  the  unperturbed  asteroid  orbit  about  the  impact  point  must  be 
mapped  onto  the  corresponding  positions  and  times  of  the  perturbed  asteroid  orbit.  That 
is,  a  one-for-one  correlation  between  where  the  asteroid  would  have  been  if  not  for  the 
perturbation  and  where  the  asteroid  is  due  to  the  perturbation  must  be  developed.  This  is 
the  key  step  in  deteirnining  the  effect  of  the  impulse  perturbation. 

The  true  anomalies  and  time  from  perihelion  of  the  perturbed  and  unperturbed 
orbits  are  related  by  Kepler's  Equation  and  not  a  simple  function.  The  approach  used  to 
achieve  the  required  mapping  is  as  follows: 

a)  Unperturbed  Orbital  Positions 

Develop  an  evenly  spaced  window,  or  mesh,  of  time  about  the  impact 
position  wide  enough  to  include  the  perturbed  orbits'  minimum  Earth  separation.  (A  first 
estimate  of  the  required  width  of  this  mesh  is  achieved  by  a  numerical  simulation.  From 
the  numerical  simulation  and  model  development  an  interval  width  of  ±  1.5  x  Av  x 
Impulse  TimeimpactTime  was  used  to  provide  a  sufficiently  wide  mesh  to  obtain  a  solution 
without  excessive  computation  time.)  By  inverting  Kepler's  Equation,  the  true  anomalies 
of  the  mesh  points  may  be  determined.  From  the  true  anomalies  of  the  mesh  points  the 
orbital  positions  may  be  found  from  Equations  (3)  and  (4). 

b)  Perturbed  Orbital  Positions 

From  the  relationship  of  the  time  of  impact  (known)  with  the  time  of 
impulse  (also  known)  the  center  of  the  mesh  may  be  determined  for  the  perturbed  orbit. 
Again,  the  true  anomalies  of  the  mesh  points  and  the  orbital  positions  for  the  perturbed 
orbit  may  be  determined  by  use  of  Kepler's  Equation  and  Equations  (3)  and  (4).  This  now 
yields  the  asteroid  orbital  positions  due  to  the  perturbation. 
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c)         Earth  Orbital  Positions 

In  the  same  fashion,  the  Earth  orbital  positions  corresponding  to  the  mesh 
points  may  be  determined.  However,  since  a  circular  Earth  orbit  was  chosen,  it  is  easier  to 
relate  the  mesh  times  to  orbital  positions  by  use  of  the  Earth's  orbital  mean  motion. 

11.        Earth  to  Asteroid  Separation  Distance 

With  the  perturbed  asteroid's  orbital  positions  and  the  corresponding  Earth  orbital 
positions  known,  it  is  a  simple  matter  to  determine  the  Earth  to  asteroid  separation 
distance  at  each  of  the  mesh  points  from: 

k2 


Ri,seP  =  ^|(xi,P  -  xi,E  )2  +  (y  i,P  -  y  i,E ) 


12.  Minimum  Separation 

From  the  above  orbital  separation  distances,  the  minimum  may  be  determined.  It  is 
useful  to  express  this  separation  in  terms  of  Earth  radii.  This  enables  a  quick  evaluation  of 
whether  a  sufficient  separation  was  achieved  to  cause  a  "miss".  As  with  the  analysis 
performed  by  Ahrens  and  Harris  (1994)  the  resultant  separation  distances  scale  linearly 
with  the  applied  impulse  magnitude. 

13.  Sample  Model  Output 

A  sample  of  the  MATLAB  model  output  is  shown  in  Figures  14  and  15.  The 
impact  scenario  is  for  the  case  of  an  impact  occurring  at  a  true  anomaly  of  30°,  an  asteroid 
orbital  eccentricity  of  54,  and  a  perturbing  impulse  of  1  m  s'1  in  the  direction  of  the  velocity 
vector  occurring  0.47  asteroid  orbits  prior  to  impact.  Figure  14  displays  the  initial 
conditions  of  the  scenario  at  the  time  of  impulse.  The  impact  position  and  the  asteroid  and 
Earth  positions  at  the  time  of  impulse  are  also  displayed.  Times  prior  to  impact  are 
displayed  on  the  asteroid  orbit  in  tenths  of  an  orbital  period. 
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Figure  14.  Initial  Conditions 


Figure  15  shows  the  effect  of  the  impulse  on  the  asteroid  near  the  impact  point. 
The  unperturbed  asteroid  position  is  shown  achieving  impact  conditions  at  zero  Earth 
radii.  The  perturbed  asteroid  trajectory  is  shown  approaching  an  impact  condition,  but 
instead  reaches  a  minimum  separation  (indicated  by  the  V)  and  then  increases  in 
separation. 
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Figure  15.  Perturbed  Asteroid  to  Earth  Separation 
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The  minimum  separation  indicated  in  Figure  15  is  the  final  output  of  the  model. 
The  figures  were  generated  from  the  model  working  variables.  This  enables  the  model  to 
be  incorporated  into  a  controlling  routine  allowing  for  repeated  simulations  that  sweep 
over  ranges  of  the  input  conditions. 
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VI.       PROGRAM  VALIDATION 

Prior  to  making  an  analysis  of  various  impact  scenarios,  the  solution  method  and 
MATLAB  model  had  to  be  validated.  This  validation  was  achieved  by  use  of  a  simple 
numerical  integration  simulation  and  by  comparison  to  approximate  analytic  solutions  of 
nearly  circular  orbits. 
A.        NUMERICAL  SIMULATION 

The  numerical  integration  simulation  was  developed  usingthe  Mathworks 
MATLAB  and  SIMULINK  numerical  processors.  The  SIMULINK  diagram  and 
MATLAB  script  files  that  support  the  SIMULINK  model  are  displayed  in  Appendix  C. 

The  two-body  equation  of  motion  integrated  by  the  model  is: 

r  =  — r. 
r^ 

Expressed  as  a  two-dimensional  system  of  linear  differential  equations  in  perifocal 
coordinate  form: 

-u  x 
v    = 

X  =  X, 

v    =  ^ 
y       r2    r' 

y  =  y. 

Numerous  orbital  scenarios  were  run  to  ensure  the  numerical  integration  was 
performing  correctly.  A  fourth  order  Runge-Kutta  integration  scheme  was  used  for  the 
simulation  acting  on  a  second  order  differential  equation.  This  yields  a  solution  that  is 
analytically  exact  and  accurate  to  the  numerical  precision  of  the  computer  microprocessor. 
For  the  present  analysis,  the  computer  microprocessor  was  a  100  MHz,  32  bit,  Intel 
Pentium.  For  the  cases  chosen  for  orbital  modeling  verification,  integrating  around  one 
orbit  resulted  in  an  ending  position  the  same  as  the  starting  position  with  a  relative  error  of 
2x1 0"16  (3x1 0"7  km  error  at  1.5xl08  km,  1  AU). 
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Once  the  numerical  integration  orbit  models  were  validated,  numerous  perturbation 
simulations  were  made  to  provide  a  reference  data  base  for  the  analytic  method.  It  was 
found  that  the  analytic  method  numerical  model  was  in  excellent  agreement  with  the 
numerical  simulations.  Differences  arose  only  from  the  difference  in  mesh  size  between 
the  two  solution  methods,  with  the  analytic  solution  having  a  finer  mesh  interval. 
B.         APPROXIMATE  ANALYTIC  SOLUTIONS 

For  circular  orbits,  Ahrens  and  Harris  (1994)  performed  an  approximate  analysis 
that  yields  expressions  for  maximum  orbital  separation  due  to  impulses  applied  both 
normal  to  and  parallel  to  the  orbital  velocity  vector.  For  the  case  where  the  impulse  is 
applied  normal  to  velocity  vector,  the  maximum  separation  is  found  half  an  orbital  period 
later  with  a  magnitude  8  max  =  2AvP  /  n .  This  case  is  shown  in  Figure  1 6. 


2AvP/ti 


Perturbed 
orbit    . 


AvP/ti 


Figure  16.  Impulse  Normal  to  Velocity  Vector 

For  the  case  where  the  impulse  is  applied  parallel  to  velocity  vector  the  separation 
a  full  orbital  period  later  is  found  to  be  5  =  3AvP  per  orbit.  That  is,  for  this  case  the 
separation  increases  by  5  for  every  orbit  after  the  single  impulse.  This  case  is  shown  in 
Figure  17. 
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Figure  17.  Impulse  Parallel  to  Velocity  Vector 

Evaluation  of  these  scenarios  by  numerical  simulation  verifies  the  above 
approximate  solution.  The  separations  between  the  perturbed  and  unperturbed  orbits  at 
the  described  locations  does  indeed  agree  very  well  with  the  approximate  analytic  solution. 
Further,  evaluation  of  the  above  scenarios  by  the  analytic  method  numerical  model  using 
an  asteroid  orbital  eccentricity  of  10"5  (nearly  circular)  and  an  impulse  of  1  m  s"1  yields 
virtually  the  same  result  as  the  approximate  analytic  solution.  Table  2  summarizes  the 
performance  of  the  three  solution  methods  for  the  circular  orbit  case. 


Impulse  Direction 

Approximate 
Analytic  Solution 

Analytic  Numerical 
Method 

Numerical 
Simulation 

Normal  to  v 

3AvP=14.8 

14.8 

14.8 

Parallel  to  v 

2AvP/tc=3.15 

3.14 

3.15 

(Separation  in  Earth  Radii) 
Table  2.  Validation  of  Solution  Method 
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VII.      ANALYSIS  AND  RESULTS 
A.         ANALYSIS 

The  preceding  analytic  numerical  method  for  solution  of  the  Earth  to  asteroid 

separation  problem  was  incorporated  as  a  function  into  a  routine  that  sweeps  over  impulse 

direction  and  time  of  impulse.  In  this  manner  hundreds  of  solutions  may  be  generated  in  a 

few  minutes.  Appendix  D  contains  numerous  analyses  for  impact  scenarios  that  may 

occur  around  the  Earth's  orbit. 

1.  Long  Time  Response  Behavior 

For  a  perturbing  impulse  time  well  prior  to  impact  the  observed  behavior  resulting 
from  the  orbital  dynamics  is  not  linear.  The  Earth-to-asteroid  separation  achieved  by  an 
impulse  is  strongly  dependent  on  the  location  of  the  impulse  on  the  orbit  as  well  as  the 
direction  of  the  impulse  with  respect  to  the  orbital  velocity.  The  typical  behavior  of  an 
impulse  is  depicted  in  Figure  1 8.  Each  point  on  this  plot  represents  the  minimum 
separation  point  depicted  on  Figure  15.  In  Figure  18  the  vertical  axis  represents  the 
minimum  Earth-to-asteroid  separation  (in  Earth  radii)  that  results  from  an  impulse  of  1  m 
s"1.  The  two  horizontal  axes  represent  the  direction  of  the  impulse  and  the  time  prior  to 
impact  of  the  impulse.  The  axis  representing  direction  of  impulse  is  measured  from  0°  to 
360°  with  respect  to  the  forward  direction  of  the  velocity  vector.  The  sense  of  the 
direction  is  that  impulse  directions  between  0°  and  1 80°  point  inward  to  the  orbit  and 
directions  from  180°  to  360°  point  outward  from  the  orbit.  The  remaining  axis  represents 
impulse  time  prior  to  impact  as  a  fraction  of  the  asteroid's  original  orbital  period  (e.g.,  t/P 
=  0.5  corresponds  to  an  impulse  one  half  orbit  prior  to  impact). 
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Figure  18.  Earth  to  Asteroid  Separation 

2.  Relative  Maxima 

Impulse  times  ranging  from  0  to  VA  asteroid  orbits  prior  to  impact  yield  two 
maxima  points  as  shown  in  Figure  1 8.  Inspection  of  these  maxima  points  reveal  that  they 
occur  at  the  time  of  perihelion  passage  for  the  asteroid  approximately  one  orbit  prior  to 
impact  for  impulse  directions  parallel  and  anti-parallel  to  the  orbital  velocity  vector.  This 
relationship  holds  true  for  all  cases  considered. 

3.  Non-perihelion  Impulse  Direction 

Closer  inspection  of  Figure  18  reveals  that  for  impulse  times  less  than  that 
corresponding  to  one  orbit  prior  to  impact,  the  maximum  separation  is  achieved  if  the 
impulse  direction  is  not  aligned  either  parallel  or  anti-parallel  to  the  velocity  vector.  This 
effect  is  better  shown  in  Figure  19,  which  is  a  contour  plot  of  Figure  1 8. 


38 


Impact  True  Anomaly  =  30  deg,  Orbital  Eccentricity  =  0.5,  dv  =  1  m/s 


90  180  270  360 

Direction  of  Impulse  wrt  V  (deg)  (ccwf ) 

Figure  19.  Contour  Plot  of  Figure  18. 


A  distinct  shift  in  the  direction  of  impulse  for  maximum  separation  as  the  time  of 
impulse  becomes  closer  to  impact  can  be  seen.  This  shift  in  direction  arises  from  the  two 
ways  in  which  to  achieve  the  desired  separation.  If  sufficient  time  prior  time  prior  to 
impact  exists,  changing  the  speed  of  the  asteroid  on  its  orbit  will  shift  the  phase  of  the 
asteroid  with  respect  to  the  Earth  and  avoid  the  impact.  If  the  time  prior  to  impact  is 
short,  changing  the  direction  of  the  asteroid  laterally  with  respect  to  its  approach  to  the 
Earth  is  necessary  to  avoid  the  impact.  For  times  prior  to  impact  between  these  the 
former  and  the  latter,  a  tradeoff  between  changing  the  orbital  speed  and  displacing  the 
asteroid  laterally  on  its  orbital  path  exists  and  is  the  cause  of  the  shift  in  the  optimal 
impulse  direction. 

4.  Periodic  Growth 

Thus  far,  only  impulse  times  within  approximately  one  asteroid  orbital  period  of 
impact  have  been  considered.  Extending  the  analysis  to  beyond  this  time  period 
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demonstrates  the  manner  in  which  the  separation  grows  as  a  function  of  impulse  time. 
Figure  20  shows  this  behavior  for  the  same  impact  scenario  as  discussed  above.  Of  note  is 
the  periodic  growth  over  time  and  the  peak  displacements  occurring  at  each  perihelion 
point.  This  particular  figure  represents  a  slice  of  Figure  18  along  the  0°  impulse  direction 
extended  from  0  to  10  orbits  prior  to  impact. 

Impact  True  Anomaly  =  30  deg,  Orbital  Eccentricity  =  0.5,  dv  =  1  m/s 
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Figure  20.  Periodic  Growth  of  Separation  Distance 

RESULTS 

The  collection  of  the  numerous  impact  scenarios  modeled  by  the  above  method 
may  be  found  in  Appendix  D.  A  study  of  these  numerous  model  results  yield  the  following 
general  conclusion. 

1.  Optimum  Impulse  Condition 

Assuming  the  optimum  impulse  condition  is  achievable  in  terms  of  a  "real"  mission 
sense  (that  is,  the  booster  technology  and  energy  delivery  mechanism  exist  for  asteroid 
mitigation),  the  optimum  impulse  point  is  located  at  the  perihelion  of  the  original  asteroid 
orbit  at  least  54  orbit  prior  to  impact.  If  time  prior  to  impact  permits,  impulse  at  perihelia 
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multiple  orbits  prior  have  an  even  greater  effect  on  separating  the  asteroid  from  the  Earth 
at  the  given  impact  time.  However,  impact  prediction  becomes  a  problem  in  such  cases  as 
the  validity  of  orbital  prediction  models  in  a  general  n-body  problem  comes  into  question. 

2.  Optimum  Impulse  Direction 

If,  due  to  the  late  time  of  detection,  it  is  not  possible  to  achieve  even  the  first 
relative  maximum  for  the  optimum  impulse  condition,  then  there  exists  an  optimum 
impulse  direction  at  the  time  of  impulse  that  maximizes  the  Earth  to  object  separation.  For 
fractions  of  an  orbit  from  0.2  <  (t/P)  <  0.9  there  is  a  shift  in  the  direction  of  impulse 
toward  the  orbital  inward  normal  (for  impulses  that  increase  the  orbital  speed)  and  toward 
the  orbital  outward  normal  (for  impulses  that  decrease  the  orbital  speed)  that  yields  a 
maximum  in  separation  achievable. 

3.  Time  of  Arrival  Consideration 

Given  the  periodic  nature  of  the  impact  problem,  there  exist  conditions  where  it 
may  be  beneficial  to  delay  a  deflecting  impulse  until  an  optimal  impulse  condition  occurs. 
For  the  scenario  corresponding  to  Figure  20,  if  time  permits  delivery  of  an  impulse  two 
and  one  half  asteroid  orbits  prior  to  impact,  it  would  prove  more  advantageous  to  delay 
the  impulse  until  only  two  orbits  prior  to  impact  in  order  to  maximize  the  effect  of  the 
impulse. 

4.  Detection  Consideration 

The  difficulty  in  realizing  the  use  of  an  optimal  impulse  condition  is  that  the 
detection  of  a  colliding  object  may  occur  too  late  to  achieve  the  most  desirable  condition. 
The  earlier  the  detection  the  better  the  chance  of  deflection  with  a  much  smaller  imparted 
energy.  Unfortunately,  the  current  search  programs  and  record  of  detection  have  been 
yielding  very  short  response  times  for  NEO's  having  very  close  approaches.  The  current 
range  of  times  for  detection  has  been  on  the  order  of  five  days  prior  to  closest  approach  to 
detection  only  after  Earth  passage. 
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VIII.    APPLICATION  TO  A  REAL  CASE 

A.         IMPULSE  ACHIEVABLE 

The  analysis  of  the  energy  coupling  between  an  explosive  yield  and  an  asteroid 
performed  by  Ahrens  and  Harris  (1994)  shows  that  it  may  be  feasible  to  deflect  a  NEO 
with  an  impulse  magnitude  from  1  cm  s"1  to  a  few  m  s'1  for  a  globally  threatening  object  of 
about  1  km  in  diameter.  This  impulse  may  be  achieved  using  one  of  several  methods. 
However,  due  to  the  relatively  short  warning  times  considered  in  this  thesis,  a  nuclear  blast 
appears  to  be  the  most  efficient  energy  delivery  mechanism.  The  blast  may  be  a  standoff 
detonation  or  surface  detonation. 

For  the  standoff  detonation,  the  required  explosive  yield,  W,  may  be  determined  by 
the  approximate  expression: 

103AvD3 
W  = 


nA       ' 

from  Ahrens  and  Harris  (1994),  which  has  been  modified  to  express  yield  in  kT  of 
equivalent  TNT.  The  impulse,  Av,  is  expressed  in  m  s"1  and  the  asteroid  diameter,  D,  is  in 
km.  The  efficiency  of  neutron  production,  n,  from  the  nuclear  blast  lies  between  0.03  and 
0.3.  The  standoff  blast  efficiency  factor,  A,  is  taken  for  an  optimum  standoff  distance  of 
0.4  asteroid  radii  with  an  efficiency  of  approximately  0.3.  An  order  of  magnitude  analysis 
of  the  above  approximation  is  presented  in  Table  3. 


Impulse  (m  s"1) 

0.1  km  diameter 

1  km  diameter 

10  km  diameter 

0.01 

0.1-1  kT 

100-1000  kT 

100-1000  MT 

0.1 

1-10  kT 

1-10  MT 

1-10  GT 

1.0 

10-100  kT 

10-100  MT 

10-100  GT 

10.0 

100-1000  kT 

100-1000  MT 

100-1000  GT 

Table  3.  Impulse  and  Diameter  v.  Standoff  Explosive  Yield 


For  the  surface  detonation,  the  required  explosive  yield,  W,  may  be  determined 
from  the  approximate  expression: 

W  =  4xlO_9AvM, 


LNEO> 
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also  from  Ahrens  and  Harris  (1994)  which  again  has  been  modified  to  express  yield  in  kT 
of  equivalent  TNT.  The  impulse,  Av,  is  expressed  inms"1  and  the  asteroid  mass,  Mneo,  is 
in  kg.  An  order  of  magnitude  analysis  of  the  above  approximation  is  presented  in  Table  4. 


Impulse  (m  s"1) 

0.1  km  diameter 

1  km  diameter 

10  km  diameter 

0.01 

60  T 

60  kT 

60  MT 

0.1 

600  T 

600  kT 

600  MT 

1.0 

6kT 

6MT 

6GT 

10.0 

60  kT 

60  MT 

60  GT 

Table  4.  Impulse  and  Diameter  v.  Surface  Explosive  Yield 

B.         TOUTATIS 

Asteroid  4179  Toutatis  will  have  multiple  close  approaches  with  the  Earth  over  the 
next  15  years.  It  is  of  interest  to  apply  the  above  analysis  and  methodology  to  Toutatis  as 
if  it  were  going  to  impact  the  Earth. 

To  perform  this  analysis  it  must  be  assumed  that  the  orbit  of  Toutatis  lies  in  the 
ecliptic  plane.  This  is  not  far  from  the  true  geometry  of  Toutatis'  orbit  where  the  orbital 
inclination  is  0.47°  out  of  the  ecliptic  plane.  From  the  catalog  of  NEA  orbital  elements 
listed  compiled  by  Tholen  (1995),  the  semi-major  axis  and  eccentricity  of  Toutatis  are 
currently  listed  as  2.5154  and  0.6361,  respectively.  Assuming  that  Toutatis  and  the  Earth 
will  collide  and  that  the  Earth  is  in  a  circular  orbit  at  1  AU  yields  an  impact  at  ±38.53° 
with  respect  to  the  perihelion  of  Toutatis.  This  is  determined  by  solving  Equation  (3)  for 
true  anomaly.  Figure  21  depicts  this  relative  orientation  of  Toutatis  with  respect  to  the 
Earth.  For  the  following  analysis,  the  +38.53°  impact  location  is  chosen  for  modeling 
purposes. 
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Impact  True  Anomaly  =  38.53  deg,  Orbital  Eccentricity  =  0.6361 

-i r 

Semi-major  axis  =  2.52  AU 


Figure  21 .  Relative  Orientation  of  Toutatis  Impact 

Modeling  the  Toutatis  collision  in  the  manner  described  above  yields  the  results 
displayed  in  Figure  22.  From  the  JPL  and  NASA  Photo  Caption  (1993)  and  Press  Release 
(1996),  the  size  of  Toutatis  is  estimated  to  be  approximately  that  of  two  attached  spheres 
having  diameters  of  about  4  km  and  2.5  km.  If  both  masses  are  combined,  the  effective 
spherical  diameter  is  approximately  4.3  km.  Using  a  mean  asteroid  density  of  3000  kg  m"3 
results  in  a  mass  for  Toutatis  near  1.25xl014  kg.  From  the  previous  impulse  analysis,  an 
explosive  yield  of  about  5  MT  is  required  in  the  case  of  a  surface  detonation  and  an 
explosive  yield  from  9  to  90  MT,  depending  on  neutron  production,  is  needed  for  a 
standoff  detonation  to  achieve  a  1  cm  s"1  change  in  orbital  speed. 

Using  the  1  cms'1  orbital  speed  change  determined  above,  a  maximum  separation 
of  1 .64  Earth  radii  is  the  result  of  an  impulse  delivered  1 .02  orbits  prior  to  impact 
(perihelion  passage  of  the  prior  orbit).  While  the  separation  is  sufficient  to  cause  the 
Toutatis  to  miss  the  Earth  in  this  scenario,  a  greater  margin  of  safety  would  be  desirable. 
A  larger  separation  may  be  achieved  by  either  delivering  an  impulse  greater  than  1  cms'1 


45 


or  delivering  the  impulse  at  an  earlier  perihelion  passage.  Recognizing  the  extremely  large 
explosive  yield  requirements  for  increased  impulse  magnitudes,  it  would  appear  preferable 
to  deliver  the  impulse  at  an  earlier  time.  This  type  of  analysis  demonstrates  the  necessity 
for  detection  of  threatening  asteroids  many  orbits  prior  to  impact. 


Impact  True  Anomaly  =  38.53  deg 
Orbital  Eccentricity  =  0.6361 
dv=0.01mfe 


360 


180 
90  -aVIc*^ 


Qffec1 


«*' 


Figure  22.  Impulse  Effects  on  Toutatis-Earth  Impact 
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IX.       FURTHER  CONSIDERATIONS 

The  above  analysis  shows  promise  as  a  tool  for  rapidly  evaluating  numerous 
scenarios  for  deflecting  an  asteroid  that  is  going  to  impact  the  Earth.  The  possibility  for 
further  investigation  utilizing  this  method  presents  itself  in  analyzing  longer  response  times 
prior  to  impact,  generalizing  the  method  to  a  three-dimensional  method  and  eccentricities 
other  than  54.  Additionally,  more  work  is  needed  in  analysis  of  very  short  response  time 
impulse  effects. 

A.  LONG  RESPONSE  TIME 

The  analysis  presented  above  and  in  Appendix  D  have  been  performed  primarily 
for  impulse  times  between  0  and  1 14  orbits  prior  to  impact.  A  few  models  were  made  for 
one  impulse  direction  at  times  ranging  from  0  to  10  orbits  prior  to  impact.  However,  this 
investigation  needs  to  be  pursued  further  in  search  of  general  trends  other  than  maximum 
separations  occurring  at  perihelion  points. 

B.  THREE-DIMENSIONAL  ANALYSIS 

The  current  model  and  method  apply  only  to  two-dimensional  scenarios.  It  is  of 
interest  and  merit  to  further  generalize  the  analysis  to  the  three-dimensional  case.  This  will 
allow  for  orbital  inclinations  out  of  the  ecliptic  plane  and  better  simulate  a  variety  of  real 
scenarios. 

C.  ECCENTRICITY  VARIATIONS 

In  the  above  analysis,  other  than  for  Toutatis,  the  orbital  eccentricities  have  been 
maintained  at  Vz.  An  investigation  into  the  effects  of  more  circular  orbits  and  highly 
eccentric  orbits  is  in  order. 

D.  SHORT  RESPONSE  TIME 

The  method  presented  is  derived  from  a  two-body  representation  of  a  more 
complicated  physical  system.  This  method  is  not  valid  for  very  short  response  times  when 
the  impacting  object  is  within  the  Earth's  sphere  of  influence.  A  further  investigation  is 
desirable  to  assess  the  effects  of  impulses  in  the  three-body  problem  that  arises  when  the 
object  is  detected  close  to  the  Earth. 
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APPENDIX  A.  DETAILED  SOLUTION  METHOD 

Given  an  Earth-asteroid  collision  with  the  following  properties: 

A  circular  Earth  orbit  with  semi-major  axis  (radius),  a  E  =  1 AU 

(orbital  eccentricity,  eE  =0); 

an  asteroid  orbit  with  orbital  eccentricity  ev ; 

the  true  anomaly  at  time  of  impact,  V^^u  =  vimpMtfE ; 


and  the  time  of  the  perturbing  impulse,  Av ,  prior  to  impact 


Av 


Find  the  minimum  separation  of  the  asteroid  from  the  Earth  in  the  vicinity  of  the  old  impact  point. 
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Nomenclature 


Variables: 

a 

semi-major  axis 

E 

eccentric  anomaly 

e 

orbital  eccentricity 

h 

specific  angular  momentum 

F 

fraction 

N 

number 

n 

orbital  mean  motion 

P 

orbital  period 

P 

orbital  parameter 

R 

asteroid  to  Earth  separation  distance,  km 

r 

orbit  to  Sun  radial  distance 

V 

velocity 

X 

x-position 

y 

y-position 

Av 

impulse  quantity 

a 

angle  w.r.t.  x-coordinate  axis 

H 

gravitational  mass  parameter 

V 

true  anomaly 

P 

asteroid  to  Earth  separation  distance,  Earth  radii 

Av 

shift  in  periapsis  direction 

fraction  of  orbit  from  periapsis 
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Subscripts: 

days 

number  of  days 

E 

Earth 

i 

i*  element 

inc 

increment 

impact 

impact  position 

min 

minimum 

orbit 

orbit  count 

P 

perturbed  asteroid  orbit 

P 

periapsis 

pos 

position 

Range 

set  of  all  i's 

sep 

separation 

Sun 

Sun  parameter 

U 

unperturbed  asteroid  orbit 

X 

x-component 

y 

y-component 

Av 

impulse 

II 

parallel  component 

X 

normal  component 

Symbols: 


A 
8 


incremental  amount 
incremental  amount 
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Elliptical  Orbit  Relations 
Mean  motion: 


-# 


Period: 


n 


Radial  distance  to  gravitational  focus  (general): 


nl. 

-ei 

r  — 

1  +  ecosv 

Periapsis  radius 

*(i 

+  ecosv1) 

r    = — — 

p 

1  +  e 

Semi-major  axis: 

a"l-e 

Eccentric  anomaly  as  a  function  of  true  anomaly: 

1 

e  +  cos  v 

1  +  ecosv 

True  anomaly  as  a  function  of  eccentric  anomaly: 

cos  E  -  e 


cosv  = 


1-ecosE 

Time  since  periapsis: 

E-esinE 
t  = 


n 
True  anomaly: 
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cosv  = 


,(l-e>) 


re  e 

Position  at  a  point  in  perifocal  coordinates: 

f  =  (r  cos  v)x  +  (r  sin  v)y 
Velocity  at  a  point  in  perifocal  coordinates: 


v=   I— [(—  sin  v)x  +  (e  +  cos  v)y  I 


Parameter  of  an  orbit: 


=  a(l-e2)  = 


h2 


Specific  angular  momentum: 

h  =  f  x  v 

Eccentricity  vector  in  perifocal  coordinates: 


1 
e  =  — 


rv 

..^ 

1  ^tl2 

P- 

V     - 

_  

V 

w 

—  r  -  (r  •  v)v 


Constants 


1AU  =  1.4959787xl08km,  astronomical  unit 


G  =  6.67259x10  20km3kg  1s  2,  gravitational  constant 


MSun  =  1.9891x10 M  kg,  mass  of  the  Sun 


HSun  =GMSun  =  1.327124399355xl0nkmV2,  Sun  gravitational  parameter 


RF  =6.37814xl03km,  Earth  radius 
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Detailed  Solution  Method 
I.  Determine  unperturbed  asteroid  orbital  elements, 
a.  Perihelion  radius: 


rF,u  = 


aE(l  +  eucosvimpactu) 


1  +  e, 


,  km 


whCTe   rimpact,U   =  aE 


b.  Semi-major  axis: 


,  km 


du 

l-e„ 

c. 

Mean  mot 
nu 

ion: 

d. 

Period: 

Pu 

2tt 
,s 

,  rad  s"1 


n, 


II.  Determine  the  conditions  at  impact  with  respect  to  the  unperturbed  asteroid  orbit, 
a.  Eccentric  anomaly  at  impact: 


E  unpack  =COS" 


1  +  eUCOSVimpact,u7 


(cos"1  principle  values  are  0  <  9  <  71 ) 
b.  Fraction  of  orbit  since  (or  prior  to)  perihelion  of  impact: 


<   t> 


vPuy 


impact,p 


"  impact  ,U        eU  Sm  ^  impact,U 

ST 


if  Vtap^u  <0, 


'0 


VluA 


<0 


impact,p 


61 


III.  Determine  the  conditions  at  impulse  with  respect  to  the  unperturbed  orbit. 
a.  Fraction  of  orbit  prior  to  (or  since)  perihelion  of  impulse: 


'O 


vP„y 


'O 


V^uA 


vP„v 


Av,p  u      impaet,p  u      Av,impact 

b.  Eccentric  anomaly  at  impulse  (solve  via  Newton-Raphson  iterative  method): 


'Av,U        *-U 


f    t   ^ 

Av,p 


Av,U 


2tt 


c.  True  anomaly  at  impulse  measured  with  respect  to  unperturbed  perifocal  coordinates: 


VAv,U  =  cos 


COsEAv,U-eU 


l-eyCOsE 


Av,U/ 


(cos"1  principle  values  are  0  <  0  <  71 ) 


PnJ 


=  N        +  F 

1N  orbit  ^  L  orbit 


Av,p 


^  - 1  <  Forbit  <  - ~  . *en  vAvU  =  2n(NoMt  - 1)  +  vAV/U 


0,27i 


if  "  2  -  Forbit  <  0 ,  then  vAV/U  =  27tNorbit  -  vAV/U 


0,2ti 


if  0  <  Forttt  <  - ,  then  vAv  u  =  2nNorbit  +  vAV/U 
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0,2ti 


if  ~  ^  Forbit  <  1 ,  then  vAvAJ  =  27t(Norbit  + 1)  -  vAv #u 


0,2tc 


d.  Distance  to  focus  (Sun)  from  impulse  location: 


km 


\1  +  eucosvAV/Uj 
e.  Velocity  components  at  impulse  with  respect  to  unperturbed  perifocal  coordinates: 


"■Sun 


x,Av,U 


<'o(l-e£) 


Sil,Viv,U.kmS"' 


v—  =  J^fcT(eu+cosv-u)'kms"' 


f.  Velocity  direction  with  respect  to  unperturbed  perifocal  coordinates 
a  =  tan_1 


y,Av,U 


\Vx/Av,U/ 


(achieved  utilizing  "atan2(y,x)"  numerical  routine) 
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IV.  Now  perturb  the  asteroid  with  respect  to  velocity  vector  direction. 

a.  Impulse  velocity  components  with  respect  to  unperturbed  perifocal  coordinates: 

Avx  =  AVj|  cos  a  -  Av±  sin  a ,  km  s"1 

Avy  =  AV||  sin  a  +  Av±  cos  a ,  km  s'1 

where,  A  v.,  is  in  the  direction  of  the  velocity  vector 

and,  Av±  is  normal  to  the  velocity  vector  in  a  right  hand  sense 

b.  Velocity  components  after  perturbation  with  respect  to  unperturbed  perifocal  coordinates: 

vx,Av/p=vx^v,u+Avx,kms-1 

Vy,Av,P   =  Vy,Av,U  +  AVy  ,  km  S"1 


Avxsina 


Av^cosa 
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V.  Determine  perturbed  asteroid  orbital  elements  from  r  and  v  . 

a.  Position  vector  with  respect  to  unperturbed  perifocal  coordinates: 

W  =  Vu  =  rX/AV/Ux  +  ry  Av  uy  +  Oz , km 

where,  rx  AvU  =  rAV/U  COS  vAvU 

md>   ry,Av,U   =  rAV/U  Sin  VAv,U 

b.  Velocity  vector  with  respect  to  unperturbed  perifocal  coordinates: 

vAv,p  =  vx,av,p*  +  vy,Av,Py  +  Oz  ,  km  s-1 

c.  Specific  angular  momentum  vector  with  respect  to  unperturbed  perifocal        coordinates: 

d.  Perturbed  orbit  parameter: 


Pp  = 


Ihp)2 

H'Sun 


=  ap(l-ep),km 


where|hp|    =  hphp 


e.  Perturbed  eccentricity  vector  with  respect  to  unperturbed  perifocal  coordinates: 


ep  = 


lSun 


H 


Sun 


Av,Pl 


rAv,ply 


rAv,P        \rAv,P  '  VAv,P/VAvJ 


where  |vAvP|2  =  vAvP  •  vAvP, 


and  IVpl=rAv,u 


f.  Perturbed  orbital  eccentricity: 


ep  ~  VeP  *eP 


g.  Perturbed  orbit  semi-major  axis: 
Pp 


aP=— — j, km 

1      ep 


h.  Perturbed  orbit  mean  motion: 
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np  = 


"Sun 


,  rad  s' 


i.  Perturbed  orbital  period: 

PP= ,s 

n„ 


VI.  Determine  the  angle  between  the  perturbed  and  unperturbed  orbital  eccentricity  vectors. 

a.  Unperturbed  orbital  eccentricity  vector: 

^u  =eux  +  0y  +  0z 

b.  Rotation  of  eccentricity  vector  due  to  impulse: 


Av  =  cos  * 


V  eFev  ) 
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VII.  Determine  the  impulse  conditions  with  respect  to  the  perturbed  orbit, 
a.  True  anomaly  at  impulse  with  respect  to  the  perturbed  orbit: 


VAv,P  =  COS 


aP(l-ep)      1 


V    rAv,pep        epy 


(cos'1  principle  values  are  0  <  0  <  7C ) 

if"1  <Forbit   <~2'then   VAv,P  =  2*(Norbit  -*)  + 


Av,P 


0,  2tc 


if~2  -  Forbit  <  0 ,  then  vAV/P  =  27iNorbit  -  vAv  p 


0,  2n 


if  0  <  ForWt  <  - ,  then  vAv  P  =  27iNorbit  +  vAV/P 


0,  2tt 


if  2  "  F°*«  <  1 ' then  Vav'p  =  27U(N°^t  +  *)  "  v 


Av,P 


71 


0,  2ti 
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b.  Eccentric  anomaly  at  impulse  with  respect  to  perturbed  orbit: 


EAv,P=COS" 


f  ep+cosvAv<p  ' 


1  +  ep  cos  v 


Av,P/ 


(cos"1  principle  values  are  0  <  0  <  7t ) 
c.  Fraction  of  orbit  of  prior  to  (or  since)  perihelion  of  impulse  with  respect  to  perturbed  orbit: 


(  ^ 


\?v) 


Av,p 


^Av,P        epSm^Av,P 
2% 


if-1<Forbit<-2'then|^" 


-(n— -i)+hr-j 


Av,p 


Av,p 


if-O^Forbit<0'then 


't^ 


\?J 


=  N       - 

1N  orbit 


't> 


Av,p 


IP, 


pS 


Av,p 


if0^ForHt<^,then 


1 

2 


't> 


vPPy 


=  N  orbit  + 


'O 


Av,p 


vi  Py 


Av,p 


1 

2 


if  o^Forbit  <l,then 


'O 


vPPy 


Av,p 


(N„.+l)-(j 


Av,p 
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VIII.  Determine  a  range  of  times,  and  hence  positions,  of  the  unperturbed  asteroid  in  the  vicinity  of  the 
impact  point  on  the  unperturbed  orbit  for  mapping  to  corresponding  times,  and  positions,  of  the  perturbed 
asteroid  on  the  perturbed  orbit. 

a.  Set  time  period: 

Ndays=±2,days 

b.  Set  number  of  positions  to  map  during  time  period: 

n^  =201 

pos 

c.  Establish  limits  of  time  interval: 


'O 


fNdavsY24hrs 


days 


Pj     I   Pv   AldayA  Ihr 


/"3600s 


d.  Number  of  increments  in  time  period: 


ninc    =npos-1 


e.  Incremental  time  step: 


r  *.  ^ 


^u- 


/  4.  "* 


=  A 


vPw 


(  1  > 


\^ne) 


f.  Now  have  n^  positions  from  -N^  to  +Ndays  in  steps  of  2Ndays/ninc  parts  of  a  day  represented 
as  fractions  of  a  complete  orbit. 


t 


> 


f 


=  -A 


Range 


uv 


< 


^t^ 


<PuJ 


<+A 


'O 


PJ 


,  in  steps  of  8 


'O 


^u/ 


g.  Center  this  range  of  positions  on  the  fraction  of  orbit  since  (or  prior  to)  perihelion  of  impact: 
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pj 


'O 


'O 


+ 


Pu>'  UV 


Range, impact  u      Range  u  '    impact,p 

h.  Eccentric  anomaly  at  each  time  coordinate  (solve  via  Newton-Raphson  iterative  method): 


't^ 


PJ 


i,Range, impact 


Ei,u-eucosEi,u 
2tc 


i.  Determine  the  true  anomalies  at  each  time  coordinate  with  respect  to  the  unperturbed  perifocal 
coordinates. 


( 


Vi,U  =  cos 


-i 


cosE^-ey 


\ 


l-eyCosE 


i,u7 


(cos1  principle  values  are  0  <  0  <  71 ) 
( 


if 


if 


if 


t 

vF„, 


Range,impact 


<  --,  then  viV  =  -2n  +  viV 


<  O 


^PuA 


<0,then  Vi„  =-viU 


i,Range,impact 


'O 


IPJ. 


>-,then  vii;  =27t-viU 


,Range/impact 

j.  Determine  the  radial  distance  to  the  focus  (Sun)  for  these  time  coordinates: 
tu(l-4) 


ri,u  = 


1  +  eyCOSViu 


,  km 


k.  Determine  the  distance  components  at  each  time  coordinate  with  respect  to  the 
unperturbed  perifocal  coordinates: 

xi,u  =riUcosviU,km 

yi/u=ri,usinviU,km 
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IX.  We  now  need  to  determine  the  true  anomaly  and  fraction  of  orbit  on  the  perturbed  orbit  that 
corresponds  to  the  true  anomaly  and  fraction  of  orbit  at  the  impact  position  on  the  unperturbed  orbit. 
Determine  the  positions  and  orbit  fractions  on  the  perturbed  orbit  that  correspond  to  the  n^ 

positions  and  orbit  fractions  of  the  unperturbed  orbit. 

a.  Fraction  of  orbit  since  (or  prior  to)  perihelion  of  the  "impact"  position  with  respect  to  the 

perturbed  orbit: 


£J    -l£J  +li4v,.P 


vPPy 


"     p/impact  r '   Av,p 

b.  Establish  the  n^  orbital  positions  and  fractions  of  orbit  on  the  perturbed  orbit  that  correspond 
exactly  to  the  n^  positions  of  the  unperturbed  orbit: 


t 


^ 


J 


< t^ 


(?„\      CO 


P      Rangcimpact 


v*-T,y 


Range 


^P 


+ 


vPpJ 


impact  ,p 


c.  Eccentric  anomaly  at  each  time  coordinate  (solve  via  Newton-Raphson  iterative  method): 


A 


J  . 


EiP  -epcosEiP 
2tc 


P      i  ,Range ,  impact 

d.  Determine  the  true  anomalies  at  each  time  coordinate  with  respect  to  the  perturbed  perifocal 
coordinates. 


vi#p  =  COS 


f  cosEiP  -ep 
l-epcosEiPy 


(cos"1  principle  values  are  0  <  0  <  71 ) 


if 


if 


—J  <--, then  ViJP  =  -2ti  +  viP 

^     i, Range,  impact 


'O 


W, 


<0,then  viP  =  -ViP 


i, Rangcimpact 
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if  4- 


\ 


>-,then  vi/P=27t-vi/P 


^     i.  Range,  impact 

e.  Shift  true  anomalies  from  perturbed  coordinates  to  unperturbed  coordinates 


Vi,P,U=Vi,P+Av 


e.  Determine  the  radial  distance  to  the  focus  (Sun)  for  these  time  coordinates: 

p(l-eP) 


1  +  eD  cos  v 


,km 


i,P 

f.  Determine  the  distance  components  at  each  time  coordinate  with  respect  to  the 
unperturbed  perifocal  coordinates: 

Xi,P  =ri/PCOSVi/P,U'km 

yi/P=riPsinvi^u,km 

X.  Now  we  need  the  positions  of  the  Earth  at  the  same  n,^  positions  surrounding  the  impact  position, 
a.  True  anomaly  of  the  Earth  at  each  of  the  n^  positions: 


'O 


Range,E  impact,E 


uv 


(Pu^e) 


Range 


b.  Determine  the  distance  components  of  the  Earth  at  each  time  coordinate  with 
to  the  unperturbed  perifocal  coordinates: 

xx,e  =aEcosviE,km 
Yi,E  =aEsinvi/E,km 


respect 
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XI.  Now  must  handle  the  special  case  of 


For  this  case,  set: 


'O 


vPu/ 


<A 


'O 


Av 


UV 


Xi,P  —  Xi,U 


yi,p  =  yi,u 


while 


^O 


^Pj 


<A 


'O 


Av 


vPT,y 


XII.  Define  the  separation  distances  of  the  perturbed  and  unperturbed  asteroid  orbital  positions  from  the 
Earth  orbital  positions. 

a.  Perturbed  asteroid  to  Earth  separation  distance: 


Rsep,p  =  J(xi,p  -xi,E)2  +(yi/P-y^)2 . 


km 


b.  Unperturbed  asteroid  to  Earth  separation  distance: 


Rsep,u  =  y(xi,u  -xi,E)2  +ijua  -Yub)2  'km 

c.  Scale  the  separation  distances  by  the  radius  of  the  Earth: 
R 


r  sep, 


'sep,P 


sep-p       R, 


,  Earth  radii 


Rsep,U 
Psep,U   =  ~^ •  Earth  radli 


73 


XIII.  Finally,  select  the  minimum  value  over  the  orbital  position  range  as  the  minimum  separation 
distance. 

a.  Minimum  perturbed  asteroid  to  Earth  separation: 

(PseP'P)min   =  mln(psep,p)'  Earth  radii 

b.  Minimum  unperturbed  asteroid  to  Earth  separation  (By  the  definition  of  the  problem  this  is 
zero,  but  this  provides  a  good  check  of  the  method.): 

(pseP,u)   .    =min(PsepU),  Earth  radi 
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APPENDIX  B.  ANALYTIC  SOLUTION  METHOD,  MATLAB  MODEL 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<>/o%%%%o/o% 

%  Arbitrary  Elliptical  Orbit  Intersecting  Circular  Earth  Orbit 

% 

%  MATLAB  Script  File  Name:  nwarborb.m 

% 

%  Author:  LT  Jeffrey  T.  Elder 

%         Naval  Postgraduate  School 

%         November  1996 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%%%%% 

%  Initialize  computational  workspace 
clear 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%%%0/o%0/o%%%%0/o%0/o%% 

%  Define  constants 

AU  =  1 .4959787e8;  %  astronomical  unit 

musun  =  1 .327124399355el  1 ;  %  gravitational  parameter  for  the  Sun 

epsilon  =  5e-6;  %  Newton-Raphson  tolerance 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%0/o%%0/o%0/o0/o%0/o%%0/o 

%  Gather  input  data  from  user 

nuimpactU  =  input('Enter  True  Anomaly  for  Impact  (deg): '); 

eU  =  input('Enter  Asteroid  eccentricity: '); 

tPUdelv  =  input('Enter  (t/P)  of  Impulse  Prior  to  Impact: '); 

delV  =  input('Enter  the  Magnitude  of  the  Impulse  dV,  (m/s): '); 

thet  =  input('Enter  the  Direction  of  the  Impulse  wrt  V,  ccw+  (deg): '); 

dwi  =  delV*cos(thet*pi/180); 

dvni  -  delV*sin(thet*pi/180); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o0/o%%%%0/o%%%% 

%  Allow  for  model  performance  evaluation 

%flops(0) 

%tic 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%0/o%%0/o0/o%0/o%% 

%  Set  Earth  orbital  elements 

aE=1.0*AU; 

nE  =  sqrt(musun/aEA3); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%% 
%  Convert  impact  true  anomaly  to  radians 
nuimpactU  =  nuimpactU*pi/180; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o°/«>0/o 

%  Determine  unperturbed  asteroid  orbital  elements 

rpU  =  aE*(H-eU*cos(nuimpactU))/(l+eU);  %  Perihelion  radius 

aU  =  rpU/(  1  -eU);  %  Semi-major  axis 

nU  =  sqrt(musun/aUA3);  %  Mean  motion 

PU  =  2*pi/nU;  %  Period 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%%0/o0/o%0/o0/o0/o0/o%0/o%%%%%0/o 

%  Determine  the  conditions  of  impact  wrt  unperturbed  orbit 

EimpactU  =...  %  ...  is  line  continuation 

acos((eU+cos(nuimpactU)). . . 

/(l+eU*cos(nuimpactU)));  %  Impact  eccentric  anomaly 

tPUpimpact  =... 

(EimpactU-eU*sin(EimpactU))/(2*pi);  %  t/P  of  impact  wrt  perihelion 

ifnuimpactU<0 

tPUpimpact  =  -tPUpimpact; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o0/o0/o%0/o%%0/o0/o%0/o%%%0/o% 

tPUpimpulse  =  tPUpimpact-tPUdelv;  %  t/P  of  impulse  wrt  perihelion 

Norbit  =  fix(tPUpimpulse);  %  Number  of  whole  orbits 

Forbit  =  tPUpimpulse  -  Norbit;  %  Fraction  of  orbits 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%%%0/o%%0/o%%0/o%0/o% 

%  Allow  for  numerical  simulation  validation 

tao  =  tPUpimpulse*PU;  %  Start  time  for  simulation 

tend  =  (tPUpimpact+0. 1 25)*PU;  %  End  time  for  simulation 


%  Guess  eccentric  anomaly 
%  allow  for  +/-  tPUpimpulse 

%  Newton-Raphson  iteration 


%  End  N-R  iteration 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o?/o%%%0/o%%%0/o%0/o0/o%%%%%%0/o0/o0/o 

%  Determine  true  anomaly  of  impulse 
EimpulseU  =  pi/4; 
tptemp  =  Forbit; 
if  tPUpimpulse  <  0 

tptemp  =  -Forbit; 
end 

true  =  0; 
while  true  —  0 

fprime  =  1  -  eU*cos(EimpulseU); 

f  =  EimpulseU  -  eU*sin(EimpulseU)  -  2*pi*tptemp; 

Elast  =  EimpulseU; 

EimpulseU  =  EimpulseU  -  f/fprime; 

if  abs(EimpulseU-Elast)  <  epsilon,  true=l;,  end 
end 

if  tPUpimpulse  <  0,  EimpulseU  =  -EimpulseU;,  end 
nuimpulseU  =... 
acos((cos(EimpulseU)-eU)./... 

(l-eU*cos(EimpulseU)));  %  True  anomaly  at  impulse 

%  If  t/P  of  impulse  exceeds 
%  one  orbit,  determine  true 
%  anomaly  for  multiple  orbits 
if  (-1  <  Forbit)  &  (Forbit  <  -0.5)  %  Conditions  to  properly 

nuimpulseU  =...  %  locate  true  anomaly 

2*pi*(Norbit- 1  )+nuimpulseU; 
elseif  (-0.5  <=  Forbit)  &  (Forbit  <  0) 

nuimpulseU  =  2*pi*Norbit  -  nuimpulseU; 
elseif  (0  <=  Forbit)  &  (Forbit  <  0.5) 

nuimpulseU  =  2*pi*Norbit  +  nuimpulseU; 
elseif  (0.5  <=  Forbit)  &  (Forbit  <  1) 

nuimpulseU  =  2*pi*(Norbit+l)  -  nuimpulseU; 
end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/o%%%o/o%o/o%%%%%o/0 

%  Determine  orbital  distances  at  time  of  impulse 
rimpulseU  =  aU*(l-eUA2)/(l+eU*cos(nuimpulseU)); 
ximpulseU  =  rimpulseU*cos(nuimpulseU); 
yimpulseU  =  rimpulseU*sin(nuimpulseU); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/o%%%%%%o/0% 

%  Perturb  the  asteroid 

pU  =  aU*(  1  -eU*eU);  %  Parameter  of  orbit 

vximpulseU=... 

-sqrt(musun/pU)*sin(nuimpulseU);  %  Unperturbed  velocities 

vyimpulseU  =... 

sqrt(musun/pU)*(eU+cos(nuimpulseU)); 
dw  =  dwi/1 000.0;  %  Impulse  parallel  to  V 

dvn  =  dvni/ 1000.0;  %  Impulse  normal  to  V 

alpha  =  atan2(vyimpulseU,vximpulseU);  %  Angle  of  V  wrt  x-axis 

vximpulseP  =... 

vximpulseU+dw*cos(alpha)-dvn*sin(alpha);  %  Perturbed  velocities 

vyimpulseP  =... 
vyimpulseU  +  dw*sin(alpha)  +  dvn*cos(alpha); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Determine  perturbed  asteroid  orbit  elements  from  R  and  V 

rimpulsePvec  =  [ximpulseU  yimpulseU  0];  %  Position  vector 

vimpulsePvec  =  [vximpulseP  vyimpulseP  0];  %  Velocity  vector 

hpvec  =  cross(rimpulsePvec, vimpulsePvec);  %  Angular  momentum  vec. 

pP  =  hpvec*hpvec7musun;  %  Parameter  of  orbit 
ePvec  =  (  (vimpulsePvec* vimpulsePvec' -musun/rimpulseU)*rimpulsePvec... 

-  (rimpulsePvec*vimpulsePvec')*vimpulsePvec  )/musun;  %  Eccentricity  vector 
eP2  =  ePvec*ePvec'; 

eP  =  sqrt(eP2);  %  Eccentricity 

aP  =  pP/(l-eP2);  %  Semi-major  axis 

nP  =  sqrt(musun/aPA3);  %  Mean  motion 

PP  =  2*pi/nP;  %  Period 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%%%%0/o%%%%%% 

%  Determine  angle  perturbed  orbit  perihelion  makes  wrt  unperturbed  orbit  perihelion 
dnu  =  acos(  (ePvec*[eU  0  0]')/(eP*eU) ); 
if  ePvec(2)  <  0,  dnu  =  -dnu;,  end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%0/o0/o%%%%0/o%0/o%0/o0/o0/o 

%  Determine  impulse  true  anomaly  wrt  perturbed  orbit 
cnuimpulseP  -  aP*(l-eP2)/(rimpulseU*eP)  -  1/eP; 

nuimpulseP  =  acos(cnuimpulseP);  %  True  anomaly  at  impulse 

%  If  t/P  of  impulse  exceeds 
%  one  orbit,  determine  true 
%  anomaly  for  multiple  orbits 
if  (-1  <  Forbit)  &  (Forbit  <  -0.5)  %  Conditions  to  properly 

nuimpulseP  =... 
2*pi*(Norbit-l)+nuimpulseP;  %  locate  true  anomaly 

elseif  (-0.5  <=  Forbit)  &  (Forbit  <  0) 

nuimpulseP  =  2*pi*Norbit  -  nuimpulseP; 
elseif  (0  <=  Forbit)  &  (Forbit  <  0.5) 
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nuimpulseP  =  2*pi*Norbit  +  nuimpulseP; 
elseif  (0.5  <=  Forbit)  &  (Forbit  <  1) 

nuimpulseP  =  2*pi*(Norbit+l)  -  nuimpulseP; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%0/o%%%0/o%0/o%0/o0/o0/o%%%0/o 

%  Determine  impulse  time  wrt  perturbed  orbit 
EimpulseP  =  acos(  (eP+cnuimpulseP)/(l+eP*cnuimpulseP)  ); 
tPPpimpulse  =... 
(EimpulseP-eP*sin(EimpulseP))/2/pi;  %  t/P  of  impulse 

%  If  t/P  of  impulse  exceeds 
%  one  orbit,  determine  t/P  of 
%  impulse  for  multiple  orbits 
if  (-1  <  Forbit)  &  (Forbit  <  -0.5)  %  Conditions  to  properly 

tPPpimpulse  =  (Norbit-l)+tPPpimpulse;  %  determine  t/P  of  impulse 

elseif  (-0.5  <=  Forbit)  &  (Forbit  <  0) 

tPPpimpulse  =  Norbit  -  tPPpimpulse; 
elseif  (0  <=  Forbit)  &  (Forbit  <  0.5) 

tPPpimpulse  =  Norbit  +  tPPpimpulse; 
elseif  (0.5  <=  Forbit)  &  (Forbit  <  1) 

tPPpimpulse  =  (Norbit+1)  -  tPPpimpulse; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%0/o%%0/o0/o% 

ndaycoef  =  1.5; 

Ndays  =  ndaycoef*tPUdelv*delV;  %  npos  &  Ndays  define  initial 

%    mesh  size 
loop  =  2;  %  Make  two  passes  to  ensure 

while(loop  >=  1)  %  'mesh' is  properly  sized 

npos  =  max([201  fix(100*Ndays+l)]);  %    to  give  accurate  minimum 

nuU=[];  %  Initialize  true  anomaly 

nuP=[];  %   matrices 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%0/o%%0/o%%0/o%0/o0/o 

%  Create  a  mesh  of  orbit  positions  about  the  impact  point  on  the 
%    unperturbed  orbit 

deltPU  =  Ndays*24*3600/PU;  %  Unperturbed  mesh  half  width 

tPURangeimpact  =...  %  Unperturbed  mesh  properly 

linspace(-deltPU,deltPU,npos)+tPUpimpact;%  located  in  t/P  space 

for  i  =  1  :npos  %  Determine  true  anomaly  of 

true  =  0;  %  unperturbed  mesh  points 

EU  =  tPURangeimpact(i)*pi; 
tptemp  =  tPURangeimpact(i); 

if  tPURangeimpact(i)  <  0,  tptemp  =  -tPURangeimpact(i);,  end 

while  true  ==  0  %  Newton-Raphson  iteration 

fprime  =  1  -  eU*cos(EU); 
f  =  EU  -  eU*sin(EU)  -  2*pi*tptemp; 
Elast  =  EU; 
EU  =  EU  -  f/fprime; 

if  abs(EU-Elast)  <  epsilon,  true  =  1 ;,  end 
end  %  End  N-R  iteration 

if  tPURangeimpact(i)  <  0,  EU  =  -EU;,  end 
nuU(i)  =... 
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acos((cos(EU)-eU)./(  1  -eU*cos(EU)));       %  True  anomalies  of 

if  tPURangeimpact(i)  <  -0.5 

nuU(i)  =  -2*pi  +  nuU(i); 
elseif  tPURangeimpact(i)  <  0 

nuU(i)  =  -nuU(i); 
elseif  tPURangeimpact(i)  >  0.5 

nuU(i)  =  2*pi  -  nuU(i); 
end 
end 

rU  -  aU*(l-eUA2) ./  (l+eU*cos(nuU)); 
xU  =  rU.*cos(nuU); 
yU  =  rU.*sin(nuU); 


%  unperturbed  mesh  points 
%  Ensure  true  anomalies  are 
%  properly  located 


%  Positions  of  unperturbed 

%  mesh  points  wrt  unperturbed 

%  coordinate  frame 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Map  the  mesh  of  orbit  positions  about  the  impact  point  onto  the  perturbed  orbit 
tPPpimpact  =  tPPpimpulse  +  tPUdelv*PU/PP; 
deltPP  -  Ndays*24*3600/PP; 
tPPRangeimpact  =... 

linspace(-deltPP,deltPP,npos)+tPPpimpact;%  located  in  t/P  space 
for  i  =  1  rnpos 

true  =  0; 

EP  =  tPPRangeimpact(i)*pi; 

tptemp  =  tPPRangeimpact(i); 

if  tPPRangeimpact(i)  <  0,  tptemp  =  -tPPRangeimpact(i);,  end 

while  true  =  0 


%  t/P  of  impact  wrt  perihelion 
%  Perturbed  mesh  half  width 
%  Perturbed  mesh  properly 

%  Determine  true  anomaly 
%  of  perturbed  mesh  points 


fprime  =  1  -  eP*cos(EP); 

f  =  EP  -  eP*sin(EP)  -  2*pi*tptemp; 

Elast  =  EP; 

EP  =  EP  -  f/fprime; 

if  abs(EP-Elast)  <  epsilon,  true  =  1;,  end 
end 

if  tPPRangeimpact(i)  <  0,  EP  =  -EP;,  end 
nuP(i)  =... 

acos((cos(EP)-eP)./(  1  -eP*cos(EP))); 

if  tPPRangeimpact(i)  <  -0.5 

nuP(i)  =  -2*pi  +  nuP(i); 
elseif  tPPRangeimpact(i)  <  0 

nuP(i)  =  -nuP(i); 
elseif  tPPRangeimpact(i)  >  0.5 

nuP(i)  =  2*pi  -  nuP(i); 
end 
end 
nuPU  -  nuP  +  dnu; 


rP  =  aP*(l-eP2)  ./  (l+eP*cos(nuP)); 
xP  =  rP.*cos(nuPU); 
yP  =  rP.*sin(nuPU); 


%  Newton-Raphson  iteration 


%  End  N-R  iteration 


%  True  anomalies  of  perturbed 
%  mesh  points 
%  Ensure  true  anomalies  are 
%  properly  located 


%  Perturbed  mesh  true 

%  anomalies  wrt  unperturbed 

%  coordinate  frame 

%  Positions  of  perturbed  mesh 

%  points  wrt  unperturbed 

%  coordinate  frame 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%%0/o%0/o%%%%0/o%%%%% 

%  Map  the  mesh  of  orbit  positions  about  the  impact  point  on  the  Earth's  orbit 

delnuE  =  Ndays*24*3600*nE;  %  Earth  mesh  half  width 
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nuRangeE  =... 
linspace(-delnuE,deliiuE,npos)+nuimpactU;  %  Earth  mesh  properly  located 

%    in  true  anomaly  space 
xE  =  AU*cos(nuRangeE);  %  Positions  of  Earth  mesh 

yE  =  AU*sin(nuRangeE);  %  points  wrt  unperturbed 

%  coordinate  frame 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o%%%%0/o0/o%0/o%%%% 

%  handle  special  case  of  tPUdelv  <  deltPU 
if  tPUdelv  <  deltPU 

[v,loc]  =  min(abs(tPPRangeimpact-tPUpimpact+tPUdelv)); 

xP(l:loc)  =  xU(l:loc); 

yP(l:loc)  =  yU(l:loc); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%0/o0/o%0/o%%0/o%%%%%%0/o% 

%  Determine  the  Earth  to  asteroid  separation  in  Earth  radii 
rhosepU  =... 

sqrt((xU-xE).A2  +  (yU-yE).A2)/6378. 14;  %  Unperturbed  orbit 

rhosepP  =... 

sqrt((xP-xE).A2  +  (yP-yE).A2)/6378. 14;  %  Perturbed  orbit 

[rhosepUmin,inu]  =  min(rhosepU);  %  Minimum  unperturbed 

%     separation  (must  be  zero) 
[rhosepPmin,inp]  =  min(rhosepP);  %  Minimum  perturbed  separation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%%%%%%0/o 

if  Ndays  ==  ndaycoef*tPUdelv*delV  %  Refine  mesh 

Ndays=  1.25*abs(tPURangeimpact(inp)*PU/3600/24  ... 
-tPURangeimpact(inu)*PU/3600/24); 
end 

loop  =  loop  -  1;  %  Allow  for  loop  exit 

end  %  End  while(test)  loop 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%% 

rhosepPmin  %  Display  minimum  separation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/O0/O%%%%%%%%%%0/0% 

%toc  %  Display  model  performance 

%flops 
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APPENDIX  C.  NUMERICAL  VALIDATION  MODEL 


SIMULINK  MODEL 
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In  order  to  run  the  SIMULINK  numerical  integration  model,  the  analytic  method 
model  must  first  be  run.  The  initial  conditions  for  the  integrators  are  specified  by  the 
working  variables  in  the  analytic  model. 
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MATLAB  SCRIPT  SUPPORTING  SIMULINK  MODEL 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o%%0/o0/o%%0/o0/o%%0/o%0/o0/o%0/o%%%0/o 

% 
% 

%  orbeom2d.m 

% 

%  Equations  of  motion  for  2d  orbit 

% 

%  Author:  LT  Jeffrey  T.  Elder 

%  Naval  Postgraduate  School 

%  November  1996 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%0/o0/o0/o%%%%0/o%% 

% 


function  xdot=orb_eom(x) 


mu=  1.327124399355ell; 
AU=1.49597870e8; 
r2  =  x(3)A2  +  x(4)A2; 
mur2  =  mu/r2; 
r  =sqrt(r2); 


%  Mass  parameter  for  Sun 

%  Astronomical  unit 

%  Radial  distance  squared 

%  Radial  distance 


xdot(l)  =  -mur2*x(3)/r; 
xdot(2)  =  -mur2*x(4)/r; 
xdot(3)  =  x(l); 
xdot(4)  =  x(2); 


%vx  1st  ODE 
%vy  1st  ODE 
%x  1st  ODE 
%y  1st  ODE 
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APPENDIX  D.  ANALYTIC  METHOD  SCENARIO  MODELS 

The  following  collection  of  figures  represents  the  analysis  performed  using  the 
numerical  routine  developed  to  perform  the  analytic  solution.  The  scenarios  are  evaluated 
for  impact  true  anomalies  from  -180°  to  +180°  referenced  to  the  asteroid  orbit  perihelion. 
The  "steps"  between  impact  scenarios  is  generally  in  30°  increments,  however,  for  clarity 
purposes,  the  step  size  is  reduced  to  10°  or  even  1°  intervals  at  times.  The  primary  30° 
increment  figures  occur  where  a  "surface  plot"  and  impact  scenario  plot  appear  together. 
At  other  increment  values,  two  surface  plots  appear  together. 

For  the  surface  plots  at  impact  scenarios  approaching  ±180°,  the  surface  appears 
rough  due  to  the  coarseness  of  the  step  size  in  impulse  time  and  impulse  direction.  At 
smaller  step  sizes  the  behavior  is  smooth  and  well  behaved. 
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